---
title: "Appendix for \"Triangulating the relationship between education and attitudes towards immigration\""
author: "Qinya Feng"
date: "Last compiled on `r format(Sys.time(), '%d %B, %Y')`"
output: 
  bookdown::pdf_document2:
    number_sections: TRUE
    toc: TRUE
    toc_depth: 2
fontsize: 12pt
indent: true
bibliography: references.bib
csl: apa.csl
header-includes:
    - \usepackage{setspace}\doublespacing
    - \usepackage[maxcitenames=2, natbib=true]{biblatex}
---

```{=tex}
\clearpage
\appendix
\renewcommand\thefigure{A.\arabic{figure}}
\setcounter{figure}{0}
\renewcommand\thetable{A.\arabic{table}} 
\setcounter{table}{0}
\renewcommand{\thesection}{\arabic{section}}
```


```{r knitr setup, include=FALSE}
# install.packages('tinytex')
# tinytex::install_tinytex()
# suppress code chunks in the Appendix pdf
knitr::opts_chunk$set(echo = FALSE, message=FALSE, warning=FALSE, fig.pos = "H")
options(tinytex.verbose = TRUE)

```

```{r setup}
# The following packages are required to be installed 
library(haven) #2.5.1
library(ggplot2) #3.4.0
library(dplyr) #1.0.10
library(tidyr) #1.3.0
library(psych) #2.2.9
library(pwr) #1.3.0
library(multilevelTools) # 0.1.1. multilevelTools provides a function iccMixed() to estimate ICCs based off of multilevel models.
library(kableExtra) #1.3.4. for nice table html formating, Vignettes for kableExtra here: https://cran.r-project.org/web/packages/kableExtra/vignettes/awesome_table_in_html.html
library(gridExtra) #2.3. for arranging plots
library(patchwork) #1.1.2. for arranging plots

rm(list = ls())
setwd("H:/Dofiles/Edu_IA")


```

```{r data}

# Read in EduIA.dta produced by Data_EduIA.do and Analysis_EduIA.do
data <- read_dta("H:/Dofiles/Edu_IA/EduIA.dta")

# Create a variable list for all main variables
mainvarlist <- c("im_avg", "im_irt", "refugee_r", "labor_immigration", "language_r", "support_culture", "education_year", "tertiary_edu", "treatment", "PGI_EA", "female", "age_2009")

# Create a variable list for all immigration attitude items
imvarlist <- c("refugee_r", "labor_immigration", "language_r", "support_culture")

```

# The SALTY survey, measurements, and descriptive statistics

\doublespacing

During 2009 to 2010, the Swedish Twin Registry (STR) administered the SALTY survey (Screening across the Life-Span Twin Cohort Study) to twins born between 1943 and 1958. The SALTY was a collaboration between researchers from various disciplines, including epidemiology, medicine, and social sciences. In addition to questions relevant to multiple fields, the survey also covered numerous items on social and political attitudes and behaviors. A total of 24,916 twins were contacted with paper questionnaire for the survey, and 11,372 respondents provided informed consent to have their data stored and analyzed [@magnusson_swedish_2013]. Apart from immigration attitude items and the municipality reform timing document collected by @holmlund_researchers_2007, all other individual-level variables in the main analysis included in the main text are obtained from Swedish register data.

## Measurements for education and immigration attitudes

This section provides more details on the measurements for education and immigration attitudes. Descriptions of other key variables (treatment indicator and PGI) are discussed in their respective sections.

### Education

\doublespacing

Information on education is obtained from LISA 2009 (Longitudinal integrated database for health insurance and labour market studies). Speciﬁcally, the "level" module in the Swedish education classiﬁcation system (SUN 2000) records the highest level of education that one obtained as of 2009. Years of education is conversed based on the highest level of education. Years of education is operationalised based on the standard classiﬁcation of education (ISCED) and in line with previous research with international comparison [@dinesen_estimating_2016]. Furthermore, the coding scheme is also consistent with some previous research on the Swedish education reform where details on codings are publicly available [@hjalmarsson_effect_2015; @holmlund_researchers_2007].

Specifically, the coding scheme is: (Old) compulsory school = 7 years; (New) compulsory school = 9 years; (Old) pre-senior secondary education = 9.5 years; High school = 10 or 11 or 12 years depending on program; Short university = 13 years; Longer university = 14 or 15 or 16 or 17 years depending on program, Short post-graduate degree = 18 years, Long post-graduate degree (PhD) = 20 years. In Study 1, it was mentioned that before the education reform, some schools already had 8 rather than 7 years of schooling. However, since school level information on their specific conditions at that time is not available, it is not feasible to distinguish individuals that had 8 years of schooling. Based on the Swedish education system, tertiary education corresponds to three years of post-secondary school education or longer, equivalent to 15 years of education.

### Immigration attitudes

Immigration attitude indices are based on four items in the SALTY survey: 1) Increase labour immigration to Sweden; 2) Introduce a language test to become a Swedish citizen; 3) Accept fewer refugees in Sweden; 4) Increase the economic support to immigrants so that they can preserve their own culture. Respondents were asked whether they think these policy proposals are "very bad" , "bad" , "not good or bad" , "good" , "very good". We coded "very bad" as 1, and "very good" as 5. The responses to the language test and refugee acceptance question are reversely coded, so that higher values represent more open and supportive views on immigration policy arrangement.

The Cronbach's $\alpha$ for the four items is about 0.73, which indicates that the reliability/internal consistency of the items used is acceptable [@tavakol_making_2011]. Item correlations are presented in Table \ref{tab:item-corr}.

```{r calculate alpha, results='hide'}

im_item<- data %>%
   select(all_of(imvarlist)) %>%
   na.omit() %>%
   alpha()

summary(im_item)

```

\singlespacing

```{r item-corr}

im_corr <- data %>%
   select(all_of(imvarlist)) %>%
   na.omit() %>%
   cor()

colnames(im_corr) <- c("Fewer refugees (r)", "Increase labor immigraion", "Language test (r)", "Economic support")
rownames(im_corr) <- c("Fewer refugees (reversed)", "Increase labor immigraion", "Language test (reversed)", "Economic support")

cor_table <- kbl(im_corr, align = "c", booktabs = TRUE, caption = "Item correlations", digits = 2) %>%
  kable_styling(font_size = 9, position = "center", latex_options = "HOLD_position", full_width = F) %>% 
  add_footnote("Notes: Fewer refugees (reversed) and Language test (reversed) are two items related to refugee acceptance and language tests that are reversely coded.", threeparttable = F)


cor_table


```

\doublespacing

Figure \ref{fig:measurements-attitude} shows the scree plot for the four items. It shows that there is a relatively large first factor or component, regardless of whether principal component or factor analysis is used. Factor loading shows that the variable on refugee acceptance has stronger correlation with the latent factor than other items.

\singlespacing

```{r measurements-attitude, fig.cap="Scree plot for immigration attitude items", out.width = "60%", out.height = "40%", fig.align='center'}

# Make the scree plot of the immigration attitude items
data %>%
   select(all_of(imvarlist)) %>%
   na.omit() %>%
   scree()

# Factor analysis
im_fa <- data %>%
  select(all_of(imvarlist)) %>%
  na.omit() %>% 
  fa(scores, # input data
              nfactors = 1, # number of factors
              rotate = "varimax", # rotation
              scores = "regression") # factor score estimation
#im_fa$loadings

# Create a data frame for factor loadings
loadings <- data.frame(
  Variables = character(),
  MR1 = numeric(),
  stringsAsFactors = FALSE
)

loadings <- rbind(loadings, data.frame(Variable = "Fewer refugees (reversed)", MR1 = 0.848))
loadings <- rbind(loadings, data.frame(Variable = "Increase labor immigraion", MR1 = 0.602))
loadings <- rbind(loadings, data.frame(Variable = "Language test (reversed)", MR1 = 0.561))
loadings <- rbind(loadings, data.frame(Variable = "Economic support", MR1 = 0.538))

# Add the SS loadings and Proportion Var
loadings <- rbind(loadings, data.frame(Variable = "SS loadings", MR1 = 1.686))
loadings <- rbind(loadings, data.frame(Variable = "Proportion Var", MR1 = 0.421))

# Print the loadings data frame
loadings_tab <- kbl(loadings, align = "l", booktabs = TRUE, caption = "Factor loadings") %>%
  kable_styling(position = "center", latex_options = "HOLD_position", full_width = F) %>% 
  add_footnote("Notes: Fewer refugees (reversed) and Language test (reversed) are reversely coded.", threeparttable = F)

loadings_tab


```

\doublespacing

Apart from averaging responses to construct one index, I use the graded response model to estimate an latent variable based index. The graded response model is an extension of the 2 parameter IRT model to ordinal items, where items vary in their difficulty and discrimination. I use `irt grm` in \textsc{STATA} to fit the model, and use the default method empirical Bayesian mean to estimate the variable. Figure \@ref(fig:irt-icc) is the item characteristic curve for each item and response.

```{r irt-icc, fig.cap="ICC graph for immigration policy attitude items", out.width = "80%", out.height = "60%", fig.align='center'}

knitr::include_graphics("Results/FigureA2.png")
```

Figure \@ref(fig:DV-density-plots) shows the distributions of both attitude indices in each sample. The general patterns are similar across samples.

```{r DV-density-plots, fig.cap = "Density plots for immigration attitudes in all samples.", out.width = "80%", out.height = "80%", fig.align='center'}


# Study 1: The MZ twin sample
data_MZ <- data %>%
  filter(sampleMZ == 1)

IA_MZ1 <- ggplot(aes(x= im_avg), data = data_MZ) +
    geom_density(alpha=0.8) +
  labs(x = "Avg index (MZ twins)")


IA_MZ2 <- ggplot(aes(x= im_irt), data = data_MZ) +
    geom_density(alpha=0.8) +
  labs(x = "IRT index (MZ twins)")

# Study 2: The education reform sample
data_reform <- data %>%
  filter(sampleReform == 1) 

IA_reform1 <- ggplot(aes(x= im_avg), data = data_reform) +
    geom_density(alpha=0.8) +
  labs(x = "Avg index (Reform)")

IA_reform2 <- ggplot(aes(x= im_irt), data = data_reform) +
    geom_density(alpha=0.8) +
  labs(x = "IRT index (Reform)")


# Study 3: The DZ twin sample
data_DZ <- data %>%
  filter(sampleDZ == 1)

IA_DZ1 <- ggplot(aes(x= im_avg), data = data_DZ) +
    geom_density(alpha=0.8) +
  labs(x = "Avg index (DZ twins)")

IA_DZ2 <- ggplot(aes(x= im_irt), data = data_DZ) +
    geom_density(alpha=0.8) +
  labs(x = "IRT index (DZ twins)")


reformDV_plot <- grid.arrange(IA_MZ1, IA_MZ2, IA_reform1, IA_reform2, IA_DZ1, IA_DZ2, nrow = 3, ncol = 2)



```

## Descriptive statistics

Table \@ref(tab:desc-MZ) - \@ref(tab:desc-DZ) give an overview of the samples and main variables. Observations with missing values on immigration attitude items are excluded. I use the \textsc{STATA} package \textsc{reghdfe} to absorb fixed effects in all studies. Since \textsc{reghdfe} drops singleton groups by default [@correia_reghdfe_2017;@RePEc:boc:bocode:s457874], the reform sample in Study 2 excludes those singleton observations that some birth year-municipality groups have. Twin samples in Study 1 (MZ twins) and 3 (DZ twins) are complete twin pairs.

\singlespacing

```{r desc-MZ}

# Study 1
desc_MZ <- data_MZ %>%
  select(all_of(mainvarlist)) %>%
  select(., -c("treatment", "PGI_EA")) %>%
  describe() %>%
  data.frame()

desc_MZ <- desc_MZ[c("n", "mean", "sd", "min", "max")]

# Rename row names
rownames(desc_MZ) <- c("Immigration policy view (Avg)", "Immigration policy view (IRT)", "Fewer refugees (reversed)", "Increase labor immigraion", "Language test (reversed)", "Economic support for immigrant culture", "Years of Education", "Tertiary Education", "Female", "Age in 2009")

# Descriptive statistics table
desc_MZ %>%
  kbl(digits = 2, align = "c", booktabs = TRUE, caption = "Descriptive statistics for the MZ twin sample") %>%
  kable_styling(position = "center", latex_options = "HOLD_position", full_width = F) %>% 
  add_footnote("Notes: The means and standard deviations for immigration attitude indices are slightly different from 0 and 1, as they are standardised based on responses from the whole SALTY sample.", threeparttable = F)

```
```{r desc-reform}

# Study 2
desc_reform <- data_reform %>%
  select(all_of(mainvarlist)) %>%
  select(., -PGI_EA) %>%
  describe() %>%
  data.frame()

desc_reform <- desc_reform[c("n", "mean", "sd", "min", "max")]

# Rename row names
variable_names <- c("Immigration policy view (Avg)", "Immigration policy view (IRT)", "Fewer refugees (reversed)", "Increase labor immigraion", "Language test (reversed)", "Economic support for immigrant culture", "Years of Education", "Tertiary Education", "Reform treatment", "Female", "Age in 2009")
rownames(desc_reform) <- variable_names

# Descriptive statistics table
desc_reform %>%
  kbl(digits = 2, align = "c", booktabs = TRUE, caption = "Descriptive statistics for the reform sample") %>%
  kable_styling(position = "center", latex_options = "HOLD_position", full_width = F) %>% 
  add_footnote("Notes: The means and standard deviations for immigration attitude indices are slightly different from 0 and 1, as they are standardised based on responses from the whole SALTY sample.", threeparttable = F)

```

```{r desc-DZ}

# Study 3
desc_DZ <- data_DZ %>%
  select(all_of(mainvarlist)) %>%
  select(., -treatment) %>%
  describe() %>%
  data.frame()

desc_DZ <- desc_DZ[c("n", "mean", "sd", "min", "max")]

# Rename row names
variable_names <- c("Immigration policy view (Avg)", "Immigration policy view (IRT)", "Fewer refugees (reversed)", "Increase labor immigraion", "Language test (reversed)", "Economic support for immigrant culture", "Years of Education", "Tertiary Education", "EA PGI", "Female", "Age in 2009")
rownames(desc_DZ) <- variable_names

# Descriptive statistics table
desc_DZ %>%
  kbl(digits = 2, align = "c", booktabs = TRUE, caption = "Descriptive statistics for the DZ twin sample") %>%
  kable_styling(position = "center", latex_options = "HOLD_position", full_width = F) %>% 
  add_footnote("Notes: The means and standard deviations for immigration attitude indices are slightly different from 0 and 1, as they are standardised based on responses from the whole SALTY sample.", threeparttable = F)

```

# Study 1: Discordant MZ twins

## Twin differences in educational attainment

\doublespacing

Study 1 uses a sample of 922 MZ twin pairs. Table \ref{tab:twin-diff-desc} shows descriptive statistics of twin differences in education. Figure \ref{fig:twin-diff-density-plot} is the density plot of twin differences in education. Figure \ref{fig:twin-diff-distribution-plot} further shows the number of pairs for each twin differences in education.

\singlespacing

```{r twin-diff-desc}

# Generate a dataset of MZ twin differences in key variables
diff_MZ <- data %>%
  select(c("im_avg", "im_irt", "education_year", "tertiary_edu", "LopNrParID", "TwinNr", "Zyg", "sampleMZ")) %>%
  filter(sampleMZ == 1) %>%
  na.omit() %>%
  group_by(LopNrParID) %>%
  summarize(diff_im_avg = abs(diff(im_avg)), #absolute value of twin differences
            diff_im_irt = abs(diff(im_irt)),
            diff_education_year = abs(diff(education_year)), 
            diff_tertiary_edu = abs(diff(tertiary_edu)),
            TwinNr = TwinNr,
            .groups = "drop")

desc_diffMZ <- diff_MZ %>%
  filter(TwinNr == 1) %>%
  select(c("diff_im_avg", "diff_im_irt", "diff_education_year", "diff_tertiary_edu")) %>%
  describe() %>%
  data.frame()

desc_diffMZ <- desc_diffMZ[c("n", "mean", "sd", "min", "max")]

# Rename row names
variable_names <- c("Immigration policy view (Avg)", "Immigration policy view (IRT)", "Years of Education", "Tertiary Education")
rownames(desc_diffMZ) <- variable_names

# Descriptive statistics table
desc_diffMZ %>%
  kbl(digits = 2, align = "c", booktabs = TRUE, caption = "Descriptive statistics MZ twin differences") %>%
  kable_styling(position = "center", latex_options = "HOLD_position", full_width = F) 

```

```{r twin-diff-density-plot, fig.cap = "Density plot of twin differences in educational attainment", out.width = "60%", out.height = "40%", fig.align='center'}

# Calculate twin differences in education and identify the lower level and higher level within a pair 
diffEdu_MZ <- data %>%
  select(c("education_year", "LopNrParID", "TwinNr", "Zyg", "sampleMZ")) %>%
  filter(sampleMZ == 1) %>% # Use the MZ twin sample as in Study 1
  na.omit() %>%
  group_by(LopNrParID) %>%
  summarize(diff_edu = abs(diff(education_year)), #absolute value of twin differences
            lower_edu = min(education_year), #lower level education in each pair, for the plot later
            twin_type = ifelse(education_year[1] > education_year[2], "high", "low"), #an indicator for twins with higher or lower educational attainment
            education_year = education_year, #essential variables for plotting later
            LopNrParID = LopNrParID,
            TwinNr = TwinNr,
            Zyg = Zyg,
            .groups = "drop")

# Create density plot of twin differences
density_twin <- ggplot(diffEdu_MZ, mapping = aes(x = diff_edu)) + 
  geom_density(fill = "gray", alpha = 0.7) +
  theme_bw() +
  scale_x_continuous(breaks = unique(diffEdu_MZ$diff_edu)) +
  xlab("Twin differences in educational attainment") +
  ylab("Density") +
  theme(axis.text = element_text(size = 10),
        axis.title = element_text(size = 12)) +
  ggtitle("Density plot of twin differences in educational attainment") 

density_twin


```

```{r twin-diff-distribution-plot, fig.cap="Distribution for MZ twin differences in educational attainment", out.width = "60%", out.height = "40%", fig.align='center'}

# Create a summary table for MZ twins with discordant education levels
# 494 pairs of twins are discordant in education
diff_MZ_sum <- diffEdu_MZ %>%
  filter(diff_edu != 0 & TwinNr == 1) %>%
  group_by(lower_edu, diff_edu) %>%
  summarize(Number_of_pairs = n())

# Create a heatmap
grid_MZ <- ggplot(diff_MZ_sum, aes(x = lower_edu, y = diff_edu, fill = Number_of_pairs, label = Number_of_pairs)) +
  geom_tile() +
  geom_text(size = 4, color = "black", show.legend = FALSE) + 
  scale_fill_gradient(low = "lightblue", high = "lightblue") +
  scale_y_continuous(breaks = unique(diff_MZ_sum$diff_edu)) +
  scale_x_continuous(breaks = unique(diff_MZ_sum$lower_edu)) +
  labs(title = "Distribution of MZ twin difference in Education",
       x = "Lower level of education in a MZ twin pair",
       y = "Difference in Education") +
  theme_bw() +
  theme(legend.position = "none") # 204 pairs have one twin finished upper secondary education (12 year) and the other twin has a higher level of education

grid_MZ

```

## Sensitivity power analysis

\doublespacing

Sensitivity power analysis is conducted under two sample sizes. One is the effective sample size based on the corresponding approximate intraclass correlation (ICC). The other one equals to the half of the total sample size. The results are presented in Figure \ref{fig:MZ-power-analysis}.

For the MZ twin sample (total N = 1844), the effective sample size is chosen as 1084, which is calculated based on an ICC of 0.7 (the ICC for years of education in the MZ twin sample). The halved sample size is 922. I use an alpha level of 0.05 with two-sided t tests to plot the relationship between varying effect sizes and statistical power. The results show that, even under the halved sample size, there is over 85% power to detect an effect size of 0.1 (beta). Therefore, the MZ twin study should have decent power to detect a "small" effect size under the conventional standard. The beta coefficients for years of education obtained in the study are about 0.11-0.12.

```{r MZ-ICC, results='hide'}

data %>%
  filter(sampleMZ == 1) %>%
  iccMixed(dv = "education_year",id = "LopNrParID")

```

```{r MZ-power-analysis, fig.cap="Sensitivity power analysis for the MZ twin Study", fig.width = 10, fig.height = 5, fig.align='center'}

# Set fixed parameters
# ICC based sample size
MZsample_size <- 1084  # Effective sample size of MZ twins based on ICC for years of education = 0.70
alpha <- 0.05       # 5% significance level

# Create a range of effect sizes 
effect_sizes <- seq(0.05, 0.5, by = 0.01)  # Adjust the range and step as needed

# Initialize an empty data frame to store results
MZ_results <- data.frame(EffectSize = numeric(0), Power = numeric(0))

# Loop through different effect sizes
for (i in 1:length(effect_sizes)) {
  effect_size <- effect_sizes[i]
  
  # Calculate power for the current effect size
  power <- pwr.t.test(d = effect_size, n = MZsample_size, sig.level = alpha, power = NULL, type = c("one.sample"), alternative = c("two.sided"))
  
  # Append the result to the data frame
  MZ_results <- rbind(MZ_results, data.frame(EffectSize = effect_size, Power = power$power))
}

# Create the plot
MZ_pw1 <- ggplot(MZ_results, aes(x = EffectSize, y = Power)) +
  geom_line() +
  geom_point() +
  labs(x = "Effect Size", y = "Power") +
  scale_x_continuous(breaks = seq(0, 0.5, by = 0.05)) +
  scale_y_continuous(breaks = seq(0, 1, by = 0.05)) +
  ggtitle("Sensitivity power analysis (N=1084)") +
  theme_bw() +
  theme(plot.title = element_text(size = 10))



# MZ twin differences sample size
MZdiffsample_size <- 922 # Sample size of MZ twin differences

# Initialize an empty data frame to store results
MZdiff_results <- data.frame(EffectSize = numeric(0), Power = numeric(0))

# Loop through different effect sizes
for (i in 1:length(effect_sizes)) {
  effect_size <- effect_sizes[i]
  
  # Calculate power for the current effect size
  power <- pwr.t.test(d = effect_size, n = MZdiffsample_size, sig.level = alpha, power = NULL, type = c("one.sample"), alternative = c("two.sided"))
  
  # Append the result to the data frame
  MZdiff_results <- rbind(MZdiff_results, data.frame(EffectSize = effect_size, Power = power$power))
}

# Create the plot
MZ_pw2 <- ggplot(MZdiff_results, aes(x = EffectSize, y = Power)) +
  geom_line() +
  geom_point() +
  labs(x = "Effect Size", y = "Power") +
  scale_x_continuous(breaks = seq(0, 0.5, by = 0.05)) +
  scale_y_continuous(breaks = seq(0, 1, by = 0.05)) +
  ggtitle("Sensitivity power analysis (N=922)") +
  theme_bw() +
  theme(plot.title = element_text(size = 10))

grid.arrange(MZ_pw1, MZ_pw2, ncol = 2)

```

## Considerations on control variables

\doublespacing

I include fields of study (highest degree) as an individual-level control, which has been overlooked in previous studies. Different fields of education have distinct values, structures, and resources that can influence attitudes towards immigrants differently. @stubager_education_2008 and @surridge_education_2016 argue that people educated in fields with a focus on the welfare of humans, such as arts and humanities, are likely to be more liberal than people who studied in fields focusing on managing people and manipulating objects, such as production and administration. Based on the social dominance theory, @sidanius_social_2003 primarily distinguish between college majors that imply sympathy toward lower-status and subordinate social groups, and majors that inherit sympathy and assistance toward groups with power such as business executives. Apart from economic and technical skills, fields that prepare students with interactive and communicative skills may lead them to be exposed to others' views, thereby socialising them to exchange and tolerate divergent standpoints (e.g., @van_de_werfhorst_sources_2004). Therefore, I expect that fields of study can vary in terms of their influence on the development of immigration attitude. Furthermore, individuals studying different fields may have varying considerations for their desired educational qualifications. Therefore, education orientation could potentially confound the relationship between educational attainment and immigration attitudes.

While there may be individual specific confounding factors, other suitable pre-treatment variables are not available or not sufficient in the cross-sectional survey data or registry data. Including variables such as cognitive performance in conscription records or grade variables at the 9th grade would significantly reduce the sample size, while conditioning on post-treatment factors can create spurious relationships and the bias on estimation can be in either direction (e.g., @acharya_explaining_2016).

Education orientation at one's highest level of education, as categorical variable, is identiﬁed by the "orientation" module in SUN 2000 (obtained from LISA2009, same as the educational attainment variable). Based on the official division (ﬁrst digits in the coding scheme), fields of study can be divided into 10 groups: General education; Education science and teacher training; Humanities and art; Social sciences, law and business administration; Natural science, mathematics and Information and Communication Technologies (ICTs); Engineering, manufacturing and construction; Agriculture and forestry, veterinary; Health and welfare; Services; Unknown.

## Full regression results

\noindent \textbf{The effects of educational attainment on attitudes towards immigration (corresponds to Figure 2 in the main text)} 

Table \ref{tab:MZ-im} presents results on the effects of education from within-family models and between family models. In the main text, the models include education orientation as a control variable; results from models excluding the variable are reported here.

In general, the education orientation variable does not change the results qualitatively. The estimates for education moderately fluctuate, and the within-family estimates for both years of schooling and tertiary education are higher in models without education orientation. We do see that several fields of study (compared with the reference group, General education) are significantly associated with immigration attitudes. However, only "Pedagogy and teaching" remain being significant in within family models.

\noindent \textbf{Comparison with existing research}

The effect size of interest for the MZ twin study can be considered in relation to the effect sizes of compulsory education reforms from existing research. The effects of one additional year of schooling on the attitude scale in the MZ twin study is about 0.04-0.05 standard deviation, which seems to be less considerable than the (pooled) reform effects. For example, using fuzzy RDD and pooling compulsory education reforms across 6 countries, Cavaille and Marshall (2019) found that one additional year of schooling decrease the anti-immigration attitude scales by about 0.2 standard deviations (they did not find significant effects for the Swedish reform, but noted that single country study could be underpowered). The authors interpreted the effect sizes as large, and suggested that the relatively uneducated reform-compliers (in the RDD) could be particularly susceptible to education's tolerance-inducing effects (p.261). D'Hombres and Nunziata (2016) also observed larger 2SLS estimates than the reduced form estimates. Similarly, they interpreted that there might be large effect of education among compliers - those with relatively low educational levels who are affected most positively by compulsory education reforms. They noted that the 2SLS estimates in their study may be capturing a local average treatment effect (LATE) that is larger than the average treatment effect (ATE) (p.212). In other words, one year of education could be more consequential for the subgroup of individuals who would not have pursued further education otherwise, than for individuals who would have obtained university education anyways.

However, the seemingly small effect size in the MZ twin study may not be necessarily against the theoretical expectation that higher levels of education would be more effective in promoting pro-immigration attitudes. The wide range of within-twin-pair variation also implies heterogeneity, compared with the local variation (compliers) leveraged in RDD or IV analysis. Therefore, other than potential influences from sample and country differences, the small effects observed in the twin study perhaps also to some extent speak to the existing interpretation that LATE is more sizable.

\setstretch{1.2}

```{=tex}
\begin{table}[H]
\centering
\begin{threeparttable}
\caption{The effects of education on immigration attitudes (MZ twins)}
\label{tab:MZ-im}
\scriptsize
\setlength{\tabcolsep}{1pt}
\begin{tabular}{lcccccccc} \hline
 & (1) & (2) & (3) & (4) & (5) & (6) & (7) & (8)  \\
\hline
 \textbf{Avg index} &  &  &  &  &  &  &  & \\
Years of Education & 0.140*** & 0.051* & 0.132*** & 0.056*** &  &  &  &  \\
 & (0.015) & (0.020) & (0.012) & (0.016) &  &  &  &  \\
Tertiary Education &  &  &  &  & 0.463*** & 0.097 & 0.591*** & 0.136 \\
 &  &  &  &  & (0.082) & (0.087) & (0.076) & (0.083) \\
Pedagogy and teaching & -0.012 & 0.138 &  &  & 0.452*** & 0.303* &  &  \\
 & (0.137) & (0.148) &  &  & (0.128) & (0.127) &  &  \\
Humanities and art & 0.170 & -0.142 &  &  & 0.589* & -0.008 &  &  \\
 & (0.227) & (0.197) &  &  & (0.238) & (0.203) &  &  \\
Social sciences Law Business Administration & -0.093 & -0.069 &  &  & 0.243* & 0.048 &  &  \\
 & (0.096) & (0.102) &  &  & (0.095) & (0.091) &  &  \\
Science Mathematics Data & -0.384 & -0.002 &  &  & 0.112 & 0.178 &  &  \\
 & (0.205) & (0.222) &  &  & (0.196) & (0.198) &  &  \\
Technology and Manufacturing & -0.309** & 0.072 &  &  & 0.004 & 0.176 &  &  \\
 & (0.103) & (0.107) &  &  & (0.103) & (0.101) &  &  \\
Agriculture Forestry Veterinary & -0.204 & 0.095 &  &  & 0.218 & 0.244 &  &  \\
 & (0.241) & (0.256) &  &  & (0.214) & (0.241) &  &  \\
Health care and Social care & -0.227* & -0.012 &  &  & 0.126 & 0.111 &  &  \\
 & (0.105) & (0.122) &  &  & (0.098) & (0.106) &  &  \\
Services & -0.139 & 0.042 &  &  & 0.175 & 0.169 &  &  \\
 & (0.125) & (0.132) &  &  & (0.121) & (0.120) &  &  \\
Unknown & 0.196 & 0.081 &  &  & 0.496 & 0.200 &  &  \\
 & (0.249) & (0.335) &  &  & (0.266) & (0.338) &  &  \\
Constant & -1.562*** & -0.655** & -1.670*** & -0.698*** & -0.206 & -0.152** & -0.162 & -0.046* \\
 & (0.289) & (0.210) & (0.298) & (0.191) & (0.283) & (0.059) & (0.263) & (0.018) \\
R-squared & 0.441 & 0.765 & 0.431 & 0.764 & 0.416 & 0.764 & 0.403 & 0.762 \\
 \hline
 
\textbf{IRT index} &  &  &  &  &  &  & \\
Years of Education & 0.153*** & 0.049* & 0.141*** & 0.054*** &  &  &  &  \\
 & (0.015) & (0.020) & (0.012) & (0.015) &  &  &  &  \\
Pedagogy and teaching & -0.079 & 0.125 &  &  & 0.432*** & 0.281* &  &  \\
 & (0.137) & (0.144) &  &  & (0.127) & (0.122) &  &  \\
Humanities and art & 0.100 & -0.101 &  &  & 0.560* & 0.026 &  &  \\
 & (0.209) & (0.181) &  &  & (0.221) & (0.182) &  &  \\
Social sciences Law Business Administration & -0.092 & -0.059 &  &  & 0.276** & 0.053 &  &  \\
 & (0.098) & (0.098) &  &  & (0.097) & (0.085) &  &  \\
Science Mathematics Data & -0.497* & -0.070 &  &  & 0.047 & 0.101 &  &  \\
 & (0.216) & (0.209) &  &  & (0.214) & (0.188) &  &  \\
Technology and Manufacturing & -0.300** & 0.048 &  &  & 0.043 & 0.147 &  &  \\
 & (0.106) & (0.105) &  &  & (0.106) & (0.102) &  &  \\
Agriculture Forestry Veterinary & -0.179 & 0.123 &  &  & 0.284 & 0.266 &  &  \\
 & (0.252) & (0.284) &  &  & (0.228) & (0.272) &  &  \\
Health care and Social care & -0.218* & 0.038 &  &  & 0.170 & 0.155 &  &  \\
 & (0.108) & (0.120) &  &  & (0.099) & (0.103) &  &  \\
Services & -0.129 & 0.097 &  &  & 0.215 & 0.218 &  &  \\
 & (0.125) & (0.126) &  &  & (0.122) & (0.115) &  &  \\
Unknown & 0.118 & 0.045 &  &  & 0.447 & 0.158 &  &  \\
 & (0.245) & (0.328) &  &  & (0.266) & (0.331) &  &  \\
Tertiary Education &  &  &  &  & 0.504*** & 0.095 & 0.620*** & 0.137 \\
 &  &  &  &  & (0.081) & (0.080) & (0.074) & (0.076) \\
Constant & -1.722*** & -0.626** & -1.777*** & -0.669*** & -0.236 & -0.148* & -0.163 & -0.038* \\
 & (0.291) & (0.202) & (0.294) & (0.181) & (0.285) & (0.058) & (0.260) & (0.017) \\
R-squared & 0.452 & 0.780 & 0.444 & 0.779 & 0.423 & 0.779 & 0.411 & 0.777 \\
 
Observations & 1,844 & 1,844 & 1,844 & 1,844 & 1,844 & 1,844 & 1,844 & 1,844 \\
Sex\#\#Age & YES & / & YES & / & YES & / & YES & / \\
Orientation & YES & YES & NO & NO & YES & YES & NO & NO \\
Twin FE & NO & YES & NO & YES & NO & YES & NO & YES \\ \hline
\end{tabular}
Notes: Standard errors are clustered at twin-pair level. Between-family models control on age dummies, sex, and their interaction terms, as well as birth municipality fixed effects. The reference group for education orientation is General education. *** p$<$0.001, ** p$<$0.01, * p$<$0.05.
\end{threeparttable}
\end{table}
```

\noindent \textbf{The effects of education by immigration attitude items (MZ twins)}

```{=tex}
\begin{table}[H]
\centering
\begin{threeparttable}
\caption{The effects of education by immigration attitude items (MZ twins)}
\label{tab:MZ-im-item}
\scriptsize
\setlength{\tabcolsep}{1pt}
\begin{tabular}{lcccccccc} \hline
& (1) & (2) & (3) & (4) & (5) & (6) & (7) & (8) \\
VARIABLES & refugee & refugee & labor & labor & culture & culture & language & language \\ \hline
 &  &  &  &  &  &  &  &  \\
Years of Education & 0.065** &  & 0.043 &  & 0.053* &  & 0.008 &  \\
 & (0.023) &  & (0.022) &  & (0.023) &  & (0.027) &  \\
Tertiary Education &  & 0.118 &  & 0.206* &  & -0.059 &  & 0.055 \\
 &  & (0.096) &  & (0.091) &  & (0.100) &  & (0.121) \\
Constant & 2.007*** & 2.644*** & 2.102*** & 2.506*** & 1.616*** & 2.155*** & 2.097*** & 2.170*** \\
 & (0.246) & (0.077) & (0.226) & (0.067) & (0.244) & (0.073) & (0.283) & (0.087) \\
 &  &  &  &  &  &  &  &  \\
Observations & 1,844 & 1,844 & 1,844 & 1,844 & 1,844 & 1,844 & 1,844 & 1,844 \\
R-squared & 0.764 & 0.763 & 0.661 & 0.661 & 0.660 & 0.659 & 0.651 & 0.651 \\
Sex\#\#Age & / & / & / & / & / & / & / & / \\
Orientation & YES & YES & YES & YES & YES & YES & YES & YES \\
 Twin FE & YES & YES & YES & YES & YES & YES & YES & YES \\ \hline
\end{tabular}
Notes: Results from within-family models are reported here. Dependent variables are on the original 1-5 scale. The responses to the language test and refugee acceptance
question are reversely coded. Standard errors are clustered at twin-pair level. *** p$<$0.001, ** p$<$0.01, * p$<$0.05.
\end{threeparttable}
\end{table}
```

# Study 2: The Swedish education reform

\doublespacing

In addition to a basic description of the Swedish education reform in the main text, here I provide more information and explanations about adjustments made.

The Swedish education reform has so far been employed to identify the causal effects of the reform and education in outcomes including educational attainment, earnings [@meghir_educational_2005], mortality [@meghir_education_2018], criminal conviction [@hjalmarsson_effect_2015], and political candidacy [@lindgren_can_2017]. In contrast to a null effect of comprehensive education reforms on educational attainment in some other European countries [@cavaille_education_2019], the Swedish education reform indeed increased years of schooling by about 0.2-0.7 years according to these research, under varied research questions, research designs, and sample restriction procedures. The estimates from this study also fall in this range.

\noindent \textbf{Reform assignment and sample restrictions}

To generate the reform treatment indicator, I rely on the municipality reform timing document coded by @holmlund_researchers_2007. This coding, based on various official documentations and archives, is used to match with the residence municipality of each individual. Following the recommendations provided (p.42-43) and previous practices [@hjalmarsson_effect_2015; @lindgren_can_2017], I exclude individuals in municipalities, or parishes within municipalities, for which the documentation does not clearly identify the adoption date of the reform. For example, the reform was gradually implemented within some municipalities and the exact timings for sub-regions (parishes) of these municipalities were therefore unclear. Some measurement error in the coding is expected, but compared with other coding methods, @holmlund_researchers_2007 found that the reliability ratio of this coding is high and suggests its high quality (p.45-50).

In addition, the oldest cohort in the SALTY sample were born in 1943 thus we are not able to include elder cohorts. In @holmlund_researchers_2007, the reform affected birth cohort range was in fact also identified as between 1943 and 1955 (p.8-9) to ensure the accuracy of actual residence municipality for individuals (earlier records were not accurate enough).

Some other adjustments on sample restrictions were made to further alleviate concerns about individual self selection into (or out of) the reform and conflated effects from other policy changes over time. First, in each municipality, I restrict the sample to only include individuals born within the range of 10 years before or after the first cohort affected by the reform. This window width is decided by balancing between obtaining an adequate sample size for estimations and sufficiently narrowing the window size to exclude potential biases from other policy changes. 10-year is used for obtaining the main reform sample for analysis. In Table \ref{tab:reformIM-multiwindow}, results using smaller window widths (9, 8, 7, 6 years) are also reported. Second, individuals who were one year older than the first affected birth cohort in each municipality were excluded. @hjalmarsson_effect_2015 found that the reform had significant effects on individuals who were born one year before the reform start years, which could be due to measurement error in reform assignment where some municipalities that partially introduced the reform to one year older cohorts are included; or grade repeators that were educated in the new system were coded as untreated, thus producing the positive effects on education in this specific birth cohort. Alternatively, these individuals might start schooling later than they were supposed to be [@lindgren_can_2017].

## Treatment distribution by birth year and reform cohort

This section presents some overview about the phased-in introduction of the education reform.

Table \@ref(tab:treatment-by-birth-cohorts) shows the distribution of treatment status in Study 2 by individual birth year.

\singlespacing

```{r treatment-by-birth-cohorts}
#for younger birth cohorts in the sample, a majority of them were exposed to the reform and therefore, were likely educated in the new compulsory schooling system, while for older cohorts, only a minority of them were exposed to the reform

# Descriptive statistics for the reform, by individual's birth year
reform_stats <- data %>%
  select(c("LopNr", "treatment", "firstcohort60", "reformKommun60", "sampleReform", "education_year", "im_avg", "BirthYear")) %>%
  filter(sampleReform == 1) %>%
  na.omit() %>% # remove observations with missing values
  group_by(BirthYear) %>%
  summarize("No. of individuals" = n(), "% treated" = sprintf("%.2f", mean(treatment)*100), "No. of treated individuals" = sum(treatment))

# create a table with footnote
reform_stats %>%
  na.omit() %>%
  kbl(align = "c", booktabs = TRUE, caption = "Treatment assignment by birth year") %>%
  kable_styling(position = "center", latex_options = "HOLD_position", full_width = F) %>%
  add_footnote("Notes: No.of individuals shows the number of individuals by each birth year cohort. % treated is the percentage, while No. of treated individuals is the number of individuals that were exposed to the education reform in the corresponding birth cohort.",  threeparttable = F)
  

```

\doublespacing

Table \@ref(tab:treatment-by-reform-cohorts) shows the distribution of reform treatment status in each municipality reform cohort. For municipalities that only started the reform later on, the proportions of treated individuals gradually decrease from mid reforming municipalities to late reforming municipalities.

\singlespacing

```{r treatment-by-reform-cohorts}


# get main descriptive statistics for the reform, by individual's birth year
reform_stats <- data %>%
  select(c("LopNr", "treatment", "firstcohort60", "reformKommun60", "sampleReform", "education_year", "im_avg", "BirthYear")) %>%
  filter(sampleReform == 1) %>%
  na.omit() %>% # remove observations with missing values
  group_by(firstcohort60) %>%
  summarize("No. of individuals" = n(), "% treated" = sprintf("%.2f", mean(treatment)*100), "No. of treated" = sum(treatment), "No. of municipalities" = n_distinct(reformKommun60)) %>%
  rename("Reform cohort" = firstcohort60)

# create a table with footnote
reform_stats %>%
  na.omit() %>%
  kbl(align = "c", booktabs = TRUE, caption = "Treatment assignment by municipality reform cohort") %>%
  kable_styling(position = "center", latex_options = "HOLD_position", full_width = F) %>%
  add_footnote("Notes: No.of individuals shows the number of individuals in each municipality reform cohort (i.e., the first affected cohort in these municipalities are the same). % treated is the percentage, while No. of treated is the number of individuals that were affected by the reform in the corresponding birth cohort. No. of municipalities is the number of municipalities in the corresponding reform cohort.",   threeparttable = F)

```

## Reform contents and result interpretation

\doublespacing

As mentioned in the main text, the Swedish education reform also had significant pedagogical content changes [@holmlund_researchers_2007; @hjalmarsson_effect_2015]. Notably, placement into differentiated academic and nonacademic tracks was postponed. Compared with systems with early tracking, less tracking could allow intergroup contacts among students from different social backgrounds and increase social trust [@osterman_can_2021], which can be considered as general positivity towards strangers and may contribute to more positive views on immigration [@herreros_social_2009]. In addition, English was introduced as a compulsory subject from 5th grade in the new comprehensive school system [@holmlund_researchers_2007]. Increased and diverse social interactions (e.g., inter-ethnic friendship) could be part of the package of education, while second language acquisition and exposure might also lead to exposure to more information and foreign culture as well as influence perceptions of intergroup relationships depending on contexts. These changes with potentially meaningful roles could all be contributing factors to the total effects of the reform on immigration attitudes.[^1]

[^1]: This is also one reason for why some studies also focus on the total effect of the education reform, rather than using the reform assignment as an instrument for years of schooling. The reform assignment could influence the outcome not only through years of schooling but also through the aforementioned changes (and these changes were not just the result of increased schooling years), thus violating the exclusion restriction assumption (see e.g., @meghir_education_2018).

However, the primary factor making any effect could be the increase in years of schooling, because only a minority of cohorts differ from the majority in terms of their exposure to English learning and delayed tracking. In terms of English learning, English was included as a compulsory subject in the curriculum and introduced in the 5th grade in the new comprehensive school. Before the reform, in general pupils went through grades 1 to 4 or 1 to 6 in a common school (folkskolan). In either 4th or 6th grade, more able students were selected (based on past performance) for the five or three/four year long junior-secondary school (realskolan) [@holmlund_researchers_2007]. At this time, English was already a part of the curriculum for realskola, but it was not a compulsory part of the folkskola curriculum. In the autumn of 1953, English was made compulsory in the old folksola [@hjalmarsson_effect_2015]. Thus, all treated and untreated individuals born in 1946 or later studied English starting in the 5th grade. In other words, only three oldest birth cohorts (born during 1943-1945) in the present study were subject to different exposure to English learning.

\par

In the new 9-year comprehensive school, pupils were kept together in common classes longer than in the earlier school system, i.e., delayed tracking. However, as a compromise between the opponents of early tracking and its advocates, tracking remained for the 9th grade until 1969. Pupils would follow either a vocational track, a general track, or a theoretical track preparing for upper-secondary school (Holmlund, 2007, p.3). Meanwhile, elective courses were available to students in the 7th and 8th grades, so that students could prepare themselves for the track that they desired to follow in grade 9 [@hjalmarsson_effect_2015]. In 1969, the 9th grade tracking was abandoned in favour of a completely comprehensive system [@holmlund_researchers_2007], which only allowed ability grouping of pupils in mathematics and English [@hjalmarsson_effect_2015]. Therefore, all but the youngest two birth cohorts in the present study (born between 1954-1955) were subject to either the early tracking that existed in the pre-reform period, or the intermediate system with parallel classes that allowed for tracking within the new comprehensive school. The process of postponing tracking is in contrast to the Finnish education reform in the 1970s where tracking was an integral part of the transformation of Finnish educational system [@pekkarinen_school_2009].

\par

In all, although I cannot completely rule out the possibility that factors other than years of schooling partially contribute to the total effects of the reform, their role should not be exaggerated either [@lindgren_can_2017; @hjalmarsson_effect_2015]. Moreover, even though other research designs such as RDD and IV can be used to identify the effects of increasing educational attainment, they identify local average treatment effects of educational attainment on reform compliers. Therefore, it is likely that the effects of the Swedish reform under the current research design could be more comparable with the effects under the discordant twin design, than designs such as RDD and IV.

## Sociopolitical characteristics of municipalities across reform cohorts {.tabset .tabset-pills}

\doublespacing

The identification of the reform effects rests on the assumption that, conditional on covariates, the reform timing was uncorrelated with unobservables that could also influence immigration attitudes. Therefore, a main concern is that municipalities were not randomly selected into each experimentation round. Although the selection of municipalities was balanced by the committee that assessed the reform, reform implementation contained voluntarism elements and certain criteria might weigh more in the decisions (see @holmlund_researchers_2007), thus leaving the reform assignment conditional on certain municipality characteristics.

There are several lines of evidence (from research using the reform time coding from @holmlund_researchers_2007) suggesting that this assumption should be valid. First, Hjalmarsson et al. (2015: p. 1305-1306) show that there was no alarming difference in pre-reform trends of years of schooling across municipalities. With more detailed analysis, they also show that individual's parental background does not predict reform treatment or have significant effects on years of schooling under appropriate specifications. Lindgren et al. (2017: p.225-227) also checked pre-reform trends of more municipality-level socioeconomic characteristics, finding that the average magnitudes of temporal change within reform cohorts were generally comparable. More nuanced balance tests also support the conclusion (in the Appendix of @lindgren_can_2017). Second, below I present visualisations that depict yearly sociopolitical characteristics of the 569 municipalities that appeared in Study 2, including voter turnout and the vote share of major active parties back then.^[The original election result data are from the Election Data Archive (Valdataarkivet) for the years 1948-1970 and the Swedish Electoral Data (Valdata) for the years 1911-1944, obtained from the Swedish National Data Service (SND). The material in the Election Data Archive (1948-1970) was originally collected by Göran Gustavsson from the Department of Sociology at Lund University, while the material in the Swedish Electoral Data (1911-1944) was collected by Sten Berglund from the Department of Political Science at Umeå University. To generate municipality-level aggregate results aligned with the 1960 census version of municipality identifiers/codes, as used in the municipality reform time records [@holmlund_researchers_2007],  a workflow previously used in @lindgren_access_2019 was adopted. This workflow processes the original data, which mainly pertains to lower administrative units (primarily at the parish level). It includes merging with the official 1960 census municipality codes and addressing changes in parish and municipality divisions and name changes over the period as well as other quality control measures (also see Lindgren et al. (2017)).] These visualisations reveal some differences among reform cohorts, but the temporal trends were similar, across pre and post reform periods (i.e., the disparities between reform cohorts were relatively consistent). Third, although I do not have information on municipality-level trends of shares of immigrants, the proportion of immigrants in the Swedish population was very low at the time (\<=5%) and immigrants at that time were mainly from Germany, other Scandinavian and Baltic countries. In all, the reform timing can be considered as plausibly exogenous. Furthermore, any time-invariant differences between early and late reformer groups will be accounted for by the reform cohort fixed effects.

The following graphs show the trends of sociopolitical characteristics (1940-1964) of municipalities for each reform cohort based on the means of municipal characteristics for each reform cohort at each election year (1940, 1944, 1948, 1952, 1956, 1958, 1960, 1964). The variables include voter population, turnout, and share of votes for major parties from the historical election results. These municipality-level statistics and their temporal trends are considered to capture important general municipal conditions. The staggered implementation period started from 1949 until 1962. In general, there is no clear sign that these socioeconomic conditions were associated with reform cohorts. One pattern that can be more clearly identified is that early reform cohorts had larger vote share for the Social Democratic party and lower voter share for the Agrarian party, than mid and late reform cohorts. However, the differences were relatively consistent over time.

In section \ref{addcontrol}, I also show that main model results with reform cohorts and birth cohorts fixed effects do not have substantial changes when controlling on these municipality level variables (including these municipality-level variables in the models with municipality fixed effects is not applicable).

```{r data preparation}

data_kommun <- read_dta("Kommun.dta") %>%
  as.data.frame()

# Extract municipality list from the reform sample
reform_kommun <- data %>%
  filter(sampleReform == 1) %>%
  select(reformKommun60) %>%
  distinct()

# Prepare a data frame by reform cohorts
data_kommun <- merge(data_kommun, reform_kommun, by = "reformKommun60") %>%
  filter(Ar >= 1940 & Ar <= 1964) %>%
  filter(firstcohort60 >= 1943 & firstcohort60 <= 1955) %>%
  mutate(Groups = case_when(firstcohort60 == 1943 ~ "Reform 1943",
                            firstcohort60 == 1944 ~ "Reform 1944",
                            firstcohort60 == 1945 ~ "Reform 1945",
                            firstcohort60 == 1946 ~ "Reform 1946",
                            firstcohort60 == 1947 ~ "Reform 1947",
                            firstcohort60 == 1948 ~ "Reform 1948",
                            firstcohort60 == 1949 ~ "Reform 1949",
                            firstcohort60 == 1950 ~ "Reform 1950",
                            firstcohort60 == 1951 ~ "Reform 1951",
                            firstcohort60 == 1952 ~ "Reform 1952",
                            firstcohort60 == 1953 ~ "Reform 1953",
                            firstcohort60 == 1954 ~ "Reform 1954",
                            firstcohort60 == 1955 ~ "Reform 1955",))
                        
                            
# Get aggregate means by reform cohorts  
data_kommun_agg <- data_kommun %>%
  group_by(Groups, Ar) %>%
  summarise(Votersize = mean(RostBer),
            Turnout = mean(Valdelt),
            M_share = mean(Andel_m),
            L_share = mean(Andel_fp),
            C_share = mean(Andel_c),
            S_share = mean(Andel_s),
            V_share = mean(Andel_v))

```

```{r voter-size, out.width = "60%", out.height = "30%"}

plot_votersize <- ggplot(data_kommun_agg, aes(x = Ar, y = Votersize, color = Groups)) +
  geom_line() +
  theme_bw() +
  labs(title = "Trends in eligible voter size",
       x = "Election year",
       y = "Average voter size") +
    scale_x_continuous(breaks = c(1940, 1944, 1948, 1952, 1956, 1958, 1960, 1964),
                     labels = c("1940", "1944", "1948", "1952", "1956", "1958", "1960", "1964"))

```

```{r voter-turnout, out.width = "60%", out.height = "30%"}

plot_turnout <- ggplot(data_kommun_agg, aes(x = Ar, y = Turnout, color = Groups)) +
  geom_line() +
  theme_bw() +
  labs(title = "Trends in turnout",
       x = "Election year",
       y = "Average turnout") +
    scale_x_continuous(breaks = c(1940, 1944, 1948, 1952, 1956, 1958, 1960, 1964),
                     labels = c("1940", "1944", "1948", "1952", "1956", "1958", "1960", "1964"))

```

```{r Conservative-share, out.width = "60%", out.height = "30%"}

plot_Mshare <- ggplot(data_kommun_agg, aes(x = Ar, y = M_share, color = Groups)) +
  geom_line() +
  theme_bw() +
  labs(title = "Trends in share of votes for the Conservative party",
       x = "Election year",
       y = "Average share of votes for the Conservative party") +
    scale_x_continuous(breaks = c(1940, 1944, 1948, 1952, 1956, 1958, 1960, 1964),
                     labels = c("1940", "1944", "1948", "1952", "1956", "1958", "1960", "1964"))

```

```{r Liberal-share, out.width = "60%", out.height = "30%"}

plot_Lshare <- ggplot(data_kommun_agg, aes(x = Ar, y = L_share, color = Groups)) +
  geom_line() +
  theme_bw() +
  labs(title = "Trends in share of votes for the Liberal party",
       x = "Election year",
       y = "Average share of votes for the Liberal party") +
    scale_x_continuous(breaks = c(1940, 1944, 1948, 1952, 1956, 1958, 1960, 1964),
                     labels = c("1940", "1944", "1948", "1952", "1956", "1958", "1960", "1964"))

```

```{r Agrarian-share, out.width = "60%", out.height = "30%"}

plot_Cshare <- ggplot(data_kommun_agg, aes(x = Ar, y = C_share, color = Groups)) +
  geom_line() +
  theme_bw() +
  labs(title = "Trends in share of votes for the Agrarian party",
       x = "Election year",
       y = "Average share of votes for the Agrarian party") +
    scale_x_continuous(breaks = c(1940, 1944, 1948, 1952, 1956, 1958, 1960, 1964),
                     labels = c("1940", "1944", "1948", "1952", "1956", "1958", "1960", "1964"))

```

```{r Social-democratis-share, out.width = "60%", out.height = "30%"}

plot_Sshare <- ggplot(data_kommun_agg, aes(x = Ar, y = S_share, color = Groups)) +
  geom_line() +
  theme_bw() +
  labs(title = "Trends in share of votes for the Social Democratic party",
       x = "Election year",
       y = "Average share of votes for the Social Democratic party") +
    scale_x_continuous(breaks = c(1940, 1944, 1948, 1952, 1956, 1958, 1960, 1964),
                     labels = c("1940", "1944", "1948", "1952", "1956", "1958", "1960", "1964"))

```

```{r Left-Communist-share, out.width = "60%", out.height = "30%"}

plot_Vshare <- ggplot(data_kommun_agg, aes(x = Ar, y = V_share, color = Groups)) +
  geom_line() +
  theme_bw() +
  labs(title = "Trends in share of votes for the Left/Communist party",
       x = "Election year",
       y = "Average share of votes for the Left/Communist party") +
    scale_x_continuous(breaks = c(1940, 1944, 1948, 1952, 1956, 1958, 1960, 1964),
                     labels = c("1940", "1944", "1948", "1952", "1956", "1958", "1960", "1964"))

```

```{r kommun-plots, fig.cap = "Sociopolitical characteristics by reform cohorts", out.width = "100%", out.height = "100%"}

(plot_votersize + plot_turnout) /
  (plot_Mshare + plot_Lshare) /
  (plot_Cshare + plot_Sshare) /
  (plot_Vshare + plot_spacer()) + plot_layout(guides = "collect") & theme(axis.title.y = element_blank()) & theme(axis.title = element_text(size = 6)) & theme(plot.title = element_text(size = 7)) & theme(axis.text = element_text(size = 6)) & theme(legend.key.size = unit(0.2, 'cm'), legend.text = element_text(size = 6), legend.title = element_text(size = 6)) 

```

## Descriptive statistics by reform treatment

This section provides descriptive statistics and Figure \@ref(fig:reform-on-education) for educational attainments by respondents' reform treatment status as shown in the main text.

\singlespacing

```{r reform-on-education, fig.cap = "Dsitribution of educational attainment by reform treatment", out.width = "60%", out.height = "60%", fig.align='center'}

# Descriptive statistics for the control group
desc_control <- data %>%
  filter(sampleReform == 1) %>%
  select(c("treatment", "education_year", "tertiary_edu", "im_avg", "im_irt", "female", "age_2009")) %>%
  na.omit() %>%
  subset(treatment == 0) %>% # subset the dataset with observations in the control group
  describe() %>%
  data.frame()

desc_control <- desc_control[c("n", "mean", "sd", "min", "max")]

# Rename row names
variable_names <- c("Reform treatment", "Years of Education", "Tertiary Education", "Immigration policy view (Avg)", "Immigration policy view (IRT)", "Female", "Age in 2009")
rownames(desc_control) <- variable_names

# Make the table
desc_control %>%
  kbl(digits = 2, align = "c", booktabs = TRUE, caption = "Descriptive statistics for the control group") %>%
  kable_styling(position = "center", latex_options = "HOLD_position", full_width = F) 


# Descriptive statistics for the treatment group
desc_treatment <- data %>%
  filter(sampleReform == 1) %>%
  select(c("treatment", "education_year", "tertiary_edu", "im_avg", "im_irt", "female", "age_2009")) %>%
  na.omit() %>%
  subset(treatment == 1) %>% # subset the dataset with observations in the treatment group
  describe() %>%
  data.frame()

desc_treatment <- desc_treatment[c("n", "mean", "sd", "min", "max")]

# Rename row names
variable_names <- c("Reform treatment", "Years of Education", "Tertiary Education", "Immigration policy view (Avg)", "Immigration policy view (IRT)", "Female", "Age in 2009")
rownames(desc_treatment) <- variable_names

# Make the table
desc_treatment %>%  
  kbl(digits = 2, align = "c", booktabs = TRUE, caption = "Descriptive statistics for the treatment group") %>%
  kable_styling(position = "center", latex_options = "HOLD_position", full_width = F) 

# Creat the plot of distribution of education by reform treatment status
dis_reform <- data %>%
  select(c("treatment", "education_year", "sampleReform")) %>%
  na.omit() %>%
  filter(sampleReform == 1) %>%
  ggplot(aes(x = factor(treatment), y = education_year)) +
  theme_bw() +
  geom_point(position = position_jitter(width = 0.3, height = 0.3), alpha = 0.3) +
  geom_boxplot(alpha = 0.2) +
  scale_y_continuous(breaks = unique(data$education_year)) +
  ylab("Years of education") +
  xlab("Control (N=2201, mean=11.54, sd=2.89)    Treated (N=1554, mean=12.35, sd=2.41) ") +
  theme(axis.ticks.x = element_blank(),
                axis.text.x=element_blank()) +
  ggtitle("Distribution of educational attainment by reform treatment")

dis_reform


```

## Reform effects on education (multiple window widths)

\doublespacing

Similar as existing research that reported the reform effects on education, the estimates in this study suggest that the reform significantly increased years of schooling, although the magnitudes of effect size vary depending on window widths for sample restriction.

It is also necessary to note that the spill-over effects of the reform on obtaining higher education are likely to be limited. Previous population-wide analysis by @holmlund_researchers_2007 suggests a small spill-over effect of the reform beyond compulsory education, lasting for up to two years after the compulsory stage but not impacting tertiary education (p.13). @cavaille_education_2019 also found no influence of education reforms on postsecondary education in their study (p.257).

\setstretch{1.2}

```{=tex}
\begin{landscape}
\begin{table}[H]
\centering
\begin{threeparttable}
\caption{Reform and years of education (mutiple window widths)}
\label{tab:reformEA-multiwindow}
\scriptsize
\setlength{\tabcolsep}{1pt}
\begin{tabular}{lccccccccccccccc}
\hline
 & (1) & (2) & (3) & (4) & (5) & (6) & (7) & (8) & (9) & (10) & (11) & (12) & (13) & (14) & (15) \\
VARIABLES & 10yr & 10yr & 10yr & 9yr & 9yr & 9yr & 8yr & 8yr & 8yr & 7yr & 7yr & 7yr & 6yr & 6yr & 6yr \\ \hline
Reform treatment & 0.732*** & 0.461** & 0.638*** & 0.733*** & 0.385* & 0.536** & 0.729*** & 0.298 & 0.491* & 0.697*** & 0.370 & 0.632** & 0.674*** & 0.196 & 0.394 \\
 & (0.134) & (0.166) & (0.188) & (0.135) & (0.170) & (0.192) & (0.135) & (0.179) & (0.206) & (0.134) & (0.190) & (0.216) & (0.134) & (0.206) & (0.239) \\
Constant & 11.451*** & 11.560*** & 11.516*** & 11.484*** & 11.629*** & 11.587*** & 11.491*** & 11.678*** & 11.621*** & 11.565*** & 11.716*** & 11.622*** & 11.618*** & 11.849*** & 11.804*** \\
 & (0.109) & (0.117) & (0.103) & (0.112) & (0.121) & (0.107) & (0.116) & (0.125) & (0.113) & (0.117) & (0.132) & (0.122) & (0.121) & (0.138) & (0.131) \\
 &  &  &  &  &  &  &  &  &  &  &  &  &  &  &  \\
Observations & 4,775 & 4,775 & 4,775 & 4,588 & 4,588 & 4,588 & 4,323 & 4,323 & 4,323 & 4,030 & 4,030 & 4,030 & 3,626 & 3,626 & 3,626 \\
R-squared & 0.029 & 0.035 & 0.270 & 0.030 & 0.035 & 0.269 & 0.029 & 0.036 & 0.273 & 0.026 & 0.031 & 0.274 & 0.027 & 0.033 & 0.287 \\
Controls & YES & YES & YES & YES & YES & YES & YES & YES & YES & YES & YES & YES & YES & YES & YES \\
Birth Cohort FE & YES & YES & YES & YES & YES & YES & YES & YES & YES & YES & YES & YES & YES & YES & YES \\
Reform Cohort FE & NO & YES & NO & NO & YES & NO & NO & YES & NO & NO & YES & NO & NO & YES & NO \\
 Municipality FE & NO & NO & YES & NO & NO & YES & NO & NO & YES & NO & NO & YES & NO & NO & YES \\ \hline
\end{tabular}
Notes: Standard errors are clustered at munitipality level. Control variables include sex. *** p$<$0.001, ** p$<$0.01, * p$<$0.05.
\end{threeparttable}
\end{table}
\end{landscape}
```


## Reform effects on immigration attitudes (multiple window widths)

\doublespacing

Table \@ref(tab:reformIM-multiwindow) presents the estimates for reform effects on immigration attitudes with multiple window widths for sample restriction. The estimates for reform treatment using a 10-year window shown in column 3 correspond to Figure 3 in the main text. As noted in the main text, the estimates for reform effects have large standard errors, and when narrowing the window width for sample restriction, the estimates from models controlling on reform cohorts or municipality fixed effects are not consistent in terms of effect direction. I consider it as mainly reflecting that the analysis is underpowered. If anything, given that a smaller window implies that the first affected cohorts get more weights, the decreasing estimates might suggest that the reform effects on the first few treated cohorts might be smaller than later treated cohorts.

\noindent \textbf{Comparison with existing research}

Two previous studies on the effects of education reform effects on immigration attitudes can be compared with Study 2 here. @cavaille_education_2019 employed the same Swedish compulsory education reform and used a fuzzy RD design to exploit the discontinuity in years of schooling when the reform started to be implemented nation-widely in a later stage. Therefore, the sample was restricted to be young enough to be affected by the country-wide coverage around 1962 (earliest cohort born in 1951), and the study identified the local average treatment effects of additional schooling at the cutoff (N range from 1561 to 4759 depending on the outcomes). The estimates for LATE were not statistically significant, but the authors also noted that country-level analysis in their study might be underpowered. @finseraas_does_2018 studied the effects of education reform in Norway where compulsory schooling was increased from 7 to 9 years and delayed tracking was also an element. Similar as the present study, they estimated the total effects of the reform within birth cohorts and municipalities. With a different composition of the attitude index and a slightly larger sample size (N = 4207), they found and argued that the reform had no causal impact on immigration attitudes given substantively small point estimates and small standard errors.

\setstretch{1.2}

```{=tex}
\begin{landscape}
\begin{table}[H]
\centering
\begin{threeparttable}
\caption{Reform and immigration attitudes (multiple window widths)}
\label{tab:reformIM-multiwindow}
\scriptsize
\begin{tabular}{lccccccccccccccc}
\hline
 & (1) & (2) & (3) & (4) & (5) & (6) & (7) & (8) & (9) & (10) & (11) & (12) & (13) & (14) & (15) \\
VARIABLES & 10yr & 10yr & 10yr & 9yr & 9yr & 9yr & 8yr & 8yr & 8yr & 7yr & 7yr & 7yr & 6yr & 6yr & 6yr \\ \hline
\textbf{Average index} &  &  &  &  &  &  &  &  &  &  &  &  &  &  &  \\
Reform treatment & 0.143** & 0.030 & 0.025 & 0.145** & -0.017 & 0.014 & 0.139** & 0.003 & 0.047 & 0.138** & 0.005 & 0.036 & 0.123* & -0.004 & -0.033 \\
 & (0.049) & (0.069) & (0.080) & (0.049) & (0.073) & (0.085) & (0.049) & (0.079) & (0.094) & (0.049) & (0.088) & (0.106) & (0.050) & (0.093) & (0.112) \\
Constant & -0.079* & -0.030 & -0.025 & -0.075* & -0.003 & -0.018 & -0.076* & -0.013 & -0.032 & -0.068 & -0.003 & -0.015 & -0.052 & 0.013 & 0.031 \\
 & (0.036) & (0.042) & (0.043) & (0.037) & (0.044) & (0.046) & (0.038) & (0.047) & (0.050) & (0.039) & (0.052) & (0.057) & (0.042) & (0.056) & (0.063) \\
 &  &  &  &  &  &  &  &  &  &  &  &  &  &  &  \\
Observations & 3,755 & 3,755 & 3,755 & 3,598 & 3,598 & 3,598 & 3,382 & 3,382 & 3,382 & 3,155 & 3,155 & 3,155 & 2,852 & 2,852 & 2,852 \\
R-squared & 0.009 & 0.015 & 0.206 & 0.010 & 0.016 & 0.208 & 0.009 & 0.016 & 0.212 & 0.009 & 0.015 & 0.223 & 0.010 & 0.018 & 0.233 \\

\hline
\textbf{IRT index} &  &  &  &  &  &  &  &  &  &  &  &  &  &  &  \\
Reform treatment & 0.136** & 0.012 & 0.016 & 0.138** & -0.035 & 0.011 & 0.132** & -0.014 & 0.048 & 0.129* & -0.009 & 0.045 & 0.115* & -0.017 & -0.012 \\
 & (0.050) & (0.071) & (0.083) & (0.050) & (0.074) & (0.088) & (0.050) & (0.081) & (0.097) & (0.050) & (0.090) & (0.108) & (0.051) & (0.094) & (0.116) \\
Constant & -0.075* & -0.021 & -0.022 & -0.070 & 0.007 & -0.015 & -0.070 & -0.002 & -0.031 & -0.058 & 0.010 & -0.014 & -0.042 & 0.027 & 0.027 \\
 & (0.036) & (0.043) & (0.044) & (0.037) & (0.045) & (0.047) & (0.038) & (0.048) & (0.052) & (0.039) & (0.053) & (0.059) & (0.042) & (0.057) & (0.065) \\
 &  &  &  &  &  &  &  &  &  &  &  &  &  &  &  \\
Observations & 3,755 & 3,755 & 3,755 & 3,598 & 3,598 & 3,598 & 3,382 & 3,382 & 3,382 & 3,155 & 3,155 & 3,155 & 2,852 & 2,852 & 2,852 \\
R-squared & 0.009 & 0.016 & 0.206 & 0.010 & 0.017 & 0.207 & 0.010 & 0.016 & 0.212 & 0.009 & 0.015 & 0.222 & 0.011 & 0.019 & 0.231 \\
\hline
Controls & YES & YES & YES & YES & YES & YES & YES & YES & YES & YES & YES & YES & YES & YES & YES \\
Birth Cohort FE & YES & YES & YES & YES & YES & YES & YES & YES & YES & YES & YES & YES & YES & YES & YES \\
Reform Cohort FE & NO & YES & NO & NO & YES & NO & NO & YES & NO & NO & YES & NO & NO & YES & NO \\
 Municipality FE & NO & NO & YES & NO & NO & YES & NO & NO & YES & NO & NO & YES & NO & NO & YES \\ \hline
\end{tabular}
Notes: Standard errors are clustered at munitipality level. Control variables include sex. *** p$<$0.001, ** p$<$0.01, * p$<$0.05.
\end{threeparttable}
\end{table}
\end{landscape}
```

## Reform effects with additional municipal controls \label{addcontrol}

\doublespacing

To evaluate whether pre-reform municipal characteristics affect the results, in Table \ref{tab:reform-kommuncontrol} I include the means and standard deviations of eligible voter population (RostBer, log transformed), turnout (Valdelt), and share of votes (Andel) for major parties (for the pre-reform period 1940-1948) as proxies for pre-reform municipal trends. For party vote share variables Andel, m denotes the conservative party, fp denotes the liberal party, c denotes the Agrarian party, s denotes the social democratic party, v denotes the left/communist party.

For immigration attitudes, only the turnout variable appear to be significantly related to immigration attitudes. Including it slightly decrease the estimates but does not substantially affect the results compared with Table \ref{tab:reformIM-multiwindow}.

\setstretch{1.2}

```{=tex}
\begin{table}[H]
\centering
\begin{threeparttable}
\caption{Reform effects with additional municipal controls}
\label{tab:reform-kommuncontrol}
\scriptsize
\begin{tabular}{lcccccc} \hline
 & (1) & (2) & (3) & (4) & (5) & (6) \\
VARIABLES & Avg index & Avg index & IRT index & IRT index & Education & Education \\ \hline
 &  &  &  &  &  &  \\
Reform treatment & 0.120* & 0.022 & 0.112* & 0.005 & 0.531*** & 0.441** \\
 & (0.049) & (0.068) & (0.050) & (0.069) & (0.130) & (0.165) \\
(mean) RostBer & -0.001 & 0.013 & -0.012 & -0.007 & 0.094 & 0.094 \\
 & (0.042) & (0.044) & (0.043) & (0.043) & (0.131) & (0.118) \\
(mean) Valdelt & -1.009 & -1.162* & -1.011* & -1.131* & -0.413 & -0.386 \\
 & (0.516) & (0.515) & (0.514) & (0.514) & (1.423) & (1.462) \\
(mean) Andel\_m & -5.115 & -6.609 & -5.372 & -6.669 & -16.914* & -17.744* \\
 & (4.169) & (4.243) & (4.058) & (4.145) & (7.469) & (8.045) \\
(mean) Andel\_fp & -5.233 & -6.650 & -5.554 & -6.781 & -18.085* & -18.601* \\
 & (4.194) & (4.219) & (4.084) & (4.119) & (7.539) & (8.066) \\
(mean) Andel\_c & -5.526 & -6.835 & -5.944 & -7.092 & -21.124** & -21.726** \\
 & (4.191) & (4.223) & (4.078) & (4.120) & (7.475) & (8.018) \\
(mean) Andel\_s & -5.501 & -6.883 & -5.910 & -7.111 & -19.167** & -19.802* \\
 & (4.148) & (4.189) & (4.037) & (4.090) & (7.389) & (7.927) \\
(mean) Andel\_v & -5.325 & -6.581 & -5.860 & -6.904 & -19.123* & -19.517* \\
 & (4.286) & (4.296) & (4.158) & (4.187) & (7.654) & (8.232) \\
(sd) RostBer & 0.670 & 0.780 & 0.628 & 0.724 & 2.135* & 2.277* \\
 & (0.411) & (0.410) & (0.396) & (0.389) & (0.941) & (0.976) \\
(sd) Valdelt & 0.165 & -0.031 & 0.097 & -0.114 & -0.784 & -1.136 \\
 & (1.168) & (1.158) & (1.138) & (1.137) & (3.158) & (3.304) \\
(sd) Andel\_m & 1.363 & 1.387 & 1.123 & 1.034 & -1.411 & -1.400 \\
 & (1.628) & (1.543) & (1.624) & (1.539) & (5.337) & (5.389) \\
(sd) Andel\_fp & -1.786 & -1.283 & -1.637 & -1.113 & -1.844 & -1.599 \\
 & (0.966) & (0.967) & (0.966) & (0.955) & (2.918) & (3.066) \\
(sd) Andel\_c & 1.194 & 0.957 & 0.920 & 0.774 & -3.009 & -2.659 \\
 & (0.950) & (0.952) & (0.933) & (0.935) & (2.977) & (2.958) \\
(sd) Andel\_s & -1.040 & -1.065 & -0.616 & -0.625 & 0.092 & -0.223 \\
 & (1.318) & (1.345) & (1.294) & (1.319) & (3.736) & (3.749) \\
(sd) Andel\_v & 2.372 & 2.110 & 2.277 & 1.874 & -1.833 & -2.783 \\
 & (1.842) & (1.834) & (1.780) & (1.757) & (4.796) & (4.834) \\
Constant & 6.036 & 7.515 & 6.456 & 7.776 & 30.563*** & 31.248*** \\
 & (4.323) & (4.367) & (4.204) & (4.262) & (7.901) & (8.478) \\
 &  &  &  &  &  &  \\
Observations & 3,755 & 3,755 & 3,755 & 3,755 & 4,775 & 4,775 \\
R-squared & 0.019 & 0.024 & 0.018 & 0.024 & 0.056 & 0.058 \\
Controls & YES & YES & YES & YES & YES & YES \\
Birth Cohort FE & YES & YES & YES & YES & YES & YES \\
Reform Cohort FE & NO & YES & NO & YES & NO & YES \\
 Municipality FE & NO & NO & NO & NO & NO & NO \\ \hline
\end{tabular}
Notes: The means and standard deviations of voter population (RostBer), turnout (Valdelt), and share of votes (Andel) for major parties (for the pre-reform period 1940-1948) are included in the model. Standard errors are clustered at munitipality level. Control variables include sex. *** p$<$0.001, ** p$<$0.01, * p$<$0.05.
\end{threeparttable}
\end{table}
```

## Reform effects by attitude items

\setstretch{1.2}

```{=tex}
\begin{table}[H]
\centering
\begin{threeparttable}
\caption{Reform effects by immigration attitude items}
\label{tab:reform-item}
\scriptsize
\setlength{\tabcolsep}{1pt}
\begin{tabular}{lcccccccccccc} \hline

 & (1) & (2) & (3) & (4) & (5) & (6) & (7) & (8) & (9) & (10) & (11) & (12) \\
VARIABLES & refugee & refugee & refugee & labor & labor & labor & culture & culture & culture & language & language & language \\ \hline
 &  &  &  &  &  &  &  &  &  &  &  &  \\
Reform treatment & 0.157** & 0.033 & 0.043 & 0.107* & -0.018 & -0.013 & 0.122** & 0.102 & 0.096 & 0.085 & -0.019 & -0.045 \\
 & (0.061) & (0.086) & (0.102) & (0.051) & (0.070) & (0.081) & (0.047) & (0.065) & (0.073) & (0.055) & (0.080) & (0.092) \\
Constant & 2.701*** & 2.756*** & 2.754*** & 2.583*** & 2.636*** & 2.634*** & 2.098*** & 2.108*** & 2.117*** & 2.335*** & 2.379*** & 2.388*** \\
 & (0.043) & (0.051) & (0.053) & (0.035) & (0.041) & (0.041) & (0.032) & (0.036) & (0.039) & (0.042) & (0.048) & (0.048) \\
 &  &  &  &  &  &  &  &  &  &  &  &  \\
Observations & 3,755 & 3,755 & 3,755 & 3,755 & 3,755 & 3,755 & 3,755 & 3,755 & 3,755 & 3,755 & 3,755 & 3,755 \\
R-squared & 0.010 & 0.016 & 0.199 & 0.007 & 0.014 & 0.199 & 0.007 & 0.010 & 0.194 & 0.011 & 0.017 & 0.188 \\
Controls & YES & YES & YES & YES & YES & YES & YES & YES & YES & YES & YES & YES \\
Birth Cohort FE & YES & YES & YES & YES & YES & YES & YES & YES & YES & YES & YES & YES \\
Reform Cohort FE & NO & YES & NO & NO & YES & NO & NO & YES & NO & NO & YES & NO \\
 Municipality FE & NO & NO & YES & NO & NO & YES & NO & NO & YES & NO & NO & YES \\ \hline
\end{tabular}
Notes: Standard errors are clustered at munitipality level. Dependent variables are on the original 1-5 scale. Control variables include sex. *** p$<$0.001, ** p$<$0.01, * p$<$0.05.
\end{threeparttable}
\end{table}
```

# Study 3: DZ twins and EA PGI

\doublespacing

In addition to a minimum explanation on PGI and the application of it in social research in the main text, here I provide a primer about GWAS, PGI, and the polygenic index repository. I also discuss sources of bias and measurement error in PGIs, the interpretations of PGIs, and some examples of using PGIs to investigate social mechanisms.

## Polygenic index

\noindent \textbf{Genome-wide association studies}

Genome-wide association study (GWAS), emerged in 2005, has allowed researchers to identify associations between thousands and millions of genetic variants throughout the genome. This means that researchers are able to actually measure and analyse genetic variants linked with traits of interest across unrelated individuals. Genetic variants here are the variations on DNA sequence between individuals. And the most common form of genetic variants in humans is single nucleotide polymorphism (SNP), where individuals differ only by one nucleotide. The nucleotide of a SNP that is more common in a population is called the major allele, and the less common one is called the minor allele. Humans receive two alleles, one from each parent. In this case, a person might have 0, 1, or 2 copies of minor allele. This number represents the genotype for a speciﬁc SNP. Genotyping thus means to identify these genetic variants that an individual possesses. Overall, humans are over 99% genetically identical in SNPs, only the remaining SNPs contains information about human variations.

In a conventional OLS-based GWAS, a speciﬁc trait of interest, such as height, is regressed onto each SNP with standard covariates like age, sex, and principal components computed from the genetic relationship matrix. This entails running millions of regression models to identify associations between SNPs and a trait. A genome-wide significant SNP needs to have a p-value (for its association with the target trait) that is lower than a multiple testing adjusted significance threshold, typically, $p < 5 \times 10^{-8}$. With developments in computing and statistical techniques and decreasing genotyping costs, GWAS research has expanded rapidly, with large sample sizes and a wide variety of traits such as educational attainment, cognitive ability, psychological traits, physical and health related traits. Large-scale GWASs often identify hundreds of or over a thousand of genetic variants that are signiﬁcantly associated with the trait in the study, with each genetic variant showing very small effect on the trait - a phenomenon called polygenicity. For example, for the ﬁrst three discovered SNPs in one of the ﬁrst GWASs on educational attainment, each explained only about 0.02% of phenotypic variance in education, which corresponds to a difference of about one month of schooling per allele [@rietveld_gwas_2013].

This theoretical and methodological development also explains why earlier studies on candidate genes often fail to be replicated. Candidate gene studies only focus on a handful of predeﬁned genes that are hypothetically relevant to the trait of interest, such as the relationship between 5HTTLPR (a genetic variant of the serotonin transporter gene) and depression [@caspi_influence_2003]. Since the inﬂuences from only a handful of genes are likely to be very small, candidate gene studies based on a few thousand individuals were dramatically underpowered to identify genetic inﬂuences. The relationship identiﬁed in those studies can thus be false positive results.

\noindent \textbf{PGI construction}

PGI represents the aggregated genetic effects across the genome by summing the alleles that an individual carries and weighting them by the allelic effect size estimates from GWASs. As a simpliﬁed example for illustration, we can consider a case of only two SNPs in the genome. As aforementioned, for each SNP a given individual can have either 0, 1 or 2 minor alleles. If we assume that GWASs on educational attainment show that each additional minor allele for the ﬁrst SNP is signiﬁcantly associated with 1 week of increase in educational attainment, and 2 weeks of increase in educational attainment for the second SNP. Then, a PGI for education can be calculated by adding the number of minor alleles for the ﬁrst SNP multiplied by 1, and the genotype for the second SNP multiplied by 2. The sum is therefore a simple measure of an individual's genetic predisposition towards the number of weeks of educational attainment.

Note, however, PGIs can be constructed in various ways although the essence of the idea remains the same. For example, PGI construction can vary in terms of the choice for GWAS discovery sample, criteria on which SNPs are included, and corrections for linkage disequilibrium (LD)\footnote{If an allele of one genetic variant is inherited or correlated with an allele of a nearby genetic variant within a given population, these alleles are in linkage disequilibrium.} (see discussions in, for example, @muslimova_rank_2023).

\noindent \textbf{Current limitations of PGI}

PGIs does not fully capture genetic variation related to target traits. The predictive power of PGI is not quite satisfactory if we compare heritability estimates from genotype data (the proportion of trait variation accounted for by measured genetic variation) with the heritability estimates from twin studies. While heritability based on all genotyped genetic variants (SNPs) has been shown to be around 20% [@rietveld_gwas_2013], a meta-analysis of twin studies reported that about 40% of the variation in educational attainment can be attributed to genetics [@branigan_variation_2013].

The main explanations for the gap between two heritability estimates are that complex traits such as social and behavioural outcomes are 1) highly polygenic (i.e., are influenced by many genetic variants), and could be contributed by 2) rare variants and 3) non-additive variance; and another explanation is that twin studies overestimate heritability (see, e.g., [@young_solving_2019]). It is expected that, as the sample sizes in GWAS increase, more comprehensive common genetic variants can be identified and incorporated in PGI construction. Genotyping rare variants and identifying non-additive variance remain challenging in the investigation of genetic architectures. However, in terms of non-additive variance, recent research has argued that genetic variance for complex traits is predominantly additive (e.g., [@hivert_estimation_2021]). Currently, the predictive power of PGIs for social and behavioural outcomes is generally lower than traits such as height and BMI. And it is acknowledged and cautioned that PGIs are poor predictors and \textit{should not be used in predicting \textbf{individual} social and behavioural outcomes} (see e.g., [@becker_resource_2021]).

\noindent \textbf{EA PGI}

Advancements in GWASs offer promise for enhancing the precision of PGIs through improvements in discovery samples and methodologies. The iterations of the polygenic index for educational attainment is one example. The polygenic index has iterated from EA1 in 2013, EA2 in 2016, EA3 in 2018, to EA4 in 2022. EA1 had a sample size of about 100,000 individuals and the polygenic index predicted about 2% of the variance in years of schooling. The most recent GWAS meta-analysis for education, EA4, has a combined sample size of about 3 million individuals and identifies 3,952 approximately uncorrelated genome-wide-significant SNPs. The EA PGI constructed explains about 12-16% of variance in education [@okbay_polygenic_2022].

The EA PGI in the polygenic index repository was constructed from a GWAS of about 1 million individuals. The incremental $R^2$ of (single-trait) EA PGI range between 6% and 11% across five validation datasets in [@becker_resource_2021] (note that the same PGI may have different predictive power in different samples if there are differences in the phenotype measure). In @dawes_polygenic_2021, the authors compared the incremental $R^2$ for EA PGI with other predictors of voter turnout in five large samples from the US and Sweden. Although the results vary across samples and different measures for turnout, the incremental $R^2$ for EA PGI in some cases are on par with variables like parental education and parental income which are often used as proxies for family socioeconomic status in social research.

## The polygenic index repository

The polygenic index repository [@becker_resource_2021] contains PGIs for 47 phenotypes and is based on meta-analysis with data from 11 participating datasets (such as the UK biobank and other registries with individuals mainly of European ancestry). It has been established with a uniform methodology and made available to help overcome obstacles faced by researcher who want to conduct GWASs and use PGIs. For example, there are administrative hurdles in accessing genetic data and summary statistics, inconsistencies of methodologies used in PGI construction, overlap between GWAS discovery sample and the target dataset).

To access the PGIs and Principal components in a dataset, researchers will need to follow the data access procedures for the specific dataset. The up-to-date procedures can be found on the Social Science Genetic Association Consortium (SSGAC) website (<https://www.thessgac.org/> pgi-repository).

The STR (genotyped individuals) is one of the participating datasets.To avoid inflation of estimates due to sample overlap, the STR cohort was excluded in the GWAS discovery stage before constructing EA PGI for individuals in the STR cohort (i.e., out-of-sample prediction). Researchers interested in using STR data must obtain approval from the Swedish Ethical Review Authority and from the Steering Committee of the Swedish Twin Registry. Researchers using STR data are required to follow the terms of a number of clauses designed to ensure protection of privacy and compliance with relevant laws. For further information please visit ki.se/en/research/the-swedish-twin-registry.

To assist interested researchers to get familiar with and use the repository, the PGI repository also provides a user guide, with more specific accounts on methodologies (and software) for PGI and principal component construction; genotyping, imputation, and phenotype definitions in Repository datasets; a list of predictive power (incremental $R^2$) of Repository PGIs in validation datasets. The authors also provide a list of suggestions on how to think through some of the interpretational issues when using PGIs.

In addition to theoretical and technical explanations, a FAQ section about the repository's research background, design, social and ethical implications, in a more accessible form, is appended in the supplemental material of @becker_resource_2021.

## Applications of PGI

Here I discuss several aspects regarding the applications and interpretations of PGI which are simplified in the main text. In the end, I also include some recent research adopting novel approaches to employ PGIs.

First, PGI is subject to measurement error and represents a noisy measure of (direct) genetic effects. The association between a genetic variant and a phenotype in most GWASs captures the effect of the variant, but likely also captures the indirect effect of the variant and confounding effects [@young_deconstructing_2019]. Indirect effects could be genetic nurture. For example, offspring genotypes can be correlated with parental genotypes (or other siblings/relatives' genotypes), while parental genotypes can influence parental behaviours which in turn exert influence on offsprings's phenotypes [@cawley_sibling_2020; @wu_estimating_2021]. Genetic nurture therefore can indirectly cause association between a variant and a phenotype. Other confounding effects could be assortative mating and population stratification. Assortative mating occurs when partners that produce offspring select each other on the basis of particular trait (e.g., education, BMI), which may induce a fraction of the association between a genetic variant and a trait to be contributed by other genetic variants that are also correlated with the trait (inflate the effect size of the genetic variant). Population stratification occurs when there are regional differences/other ancestral differences in allele frequency related to a trait of interest. If not adequately controlled for by principal components, residual population stratification can be correlated with both phenotypes and genotypes. Given these issues, family-based GWASs and other methodologies have been developed or conducted to gradually eliminate these sources of bias in GWASs [@howe_within-sibship_2022; @young_relatedness_2018]. In all, PGIs constructed based on effect sizes estimated in most GWASs captures some true genetic effects but also indirect genetic effects and other residual confounding.

Second, when applying PGIs to downstream analysis, the association between a PGI and a trait might be confounded by environmental mechanisms as discussed in the previous paragraph. But in within-family analysis, as noted, genetic inheritance from a pair of biological parents is random. Genetic differences between siblings therefore are uncorrelated with familial genetic influences that shape environments, or population stratification and assortative mating. Consequently, within-family analysis identifies the causal effects of a PGI. Since PGIs contain measurement error, within-family estimates of PGI as genetic effects are likely subject to attenuation bias [@trejo_genetic_2018; @fletcher_interpreting_2021; @becker_resource_2021].

Third, using PGIs in social research is not aimed for prediction in the same sense as some health and disease research. Despite the increasing precision of PGIs, for complex behavioural traits such as education, they tend to explain limited proportion of variance in the target traits (compared with PGIs for height and BMI, for instance). Rather, genetically informed social research has been combining PGIs with various designs (family based or unrelated individuals) to help elucidate social mechanisms, such as genetic nurture, assortative mating, and gene-environment correlations which are essentially consequential environmental and socially based mechanisms or contexts. For example, using PGI for education of a large sample of parent-offspring trios in Norway, @nivard_neither_2022 find that indirect genetic effects on children's educational attainment are primarily explained by components derived from multi-generational process of social stratification. Using within-family analysis, @fletcher_production_2023 show that the association between EA PGI and education is relatively weaker for the sibling that has higher PGI (compared with the other ones), which might indicate that compensatory processes are present for children with lower PGI. As one example for gene-environment correlation (rGE), @abdellaoui_genetic_2019 identify genetically related geographical mobility that suggest rGE. They find that individuals who are genetically predisposed to higher educational attainment in Great Britain were more likely to migrate out of lower SES areas such as coal mining areas.

## Sensitivity power analysis

Similar as the analysis for MZ twins, for the DZ twin sample (total N = 3168), the effective sample size is chosen as 2030, which is calculated based on an ICC of 0.56 (the ICC for EA PGI in the DZ twin sample, which is very close to the twin-pair ICC for education related polygenic index reported in @selzam_comparing_2019). The halved sample size is 1584. I use an alpha level of 0.05 with two-sided t tests to plot the relationship between varying effect sizes and statistical power.

The results in Figure \ref{fig:DZ-power-analysis} show that, even under the halved sample size, there is over 95% power to detect an effect size of 0.1 (beta). Therefore, the DZ twin sample should have fairly decent power to detect a small effect size. In the present DZ twin study, one sd of EA PGI corresponds to an increase of immigration attitudes scale by about 0.13-0.14 sd (equivalent to beta). So far, there is little research examining the effects of EA PGI on political attitudes, but we can compare these effect sizes with @dawes_polygenic_2021 which examines the effects of EA PGI on political participation. @dawes_polygenic_2021 shows that one sd of EA PGI increases individual average turnout across national elections by about 0.01-0.02 (with various samples), and increases the self-reported voting scale by about 0.15-0.18 sd.

```{r DZ-ICC, results='hide'}

data %>%
  filter(sampleDZ == 1) %>%
  iccMixed(dv = "PGI_EA",id = "LopNrParID")

```

```{r DZ-power-analysis, fig.cap="Sensitivity power analysis for the DZ twin Study", fig.width = 10, fig.height = 5, fig.align='center'}

# Set fixed parameters
DZsample_size <- 2030  # Effective sample size of DZ twins based on ICC of EA PGI = 0.56
alpha <- 0.05       # 5% significance level

# Create a range of effect sizes
effect_sizes <- seq(0.05, 0.5, by = 0.01)  # Adjust the range and step as needed

# Initialize an empty data frame to store results
DZresults <- data.frame(EffectSize = numeric(0), Power = numeric(0))

# Loop through different effect sizes
for (i in 1:length(effect_sizes)) {
  effect_size <- effect_sizes[i]
  
  # Calculate power for the current effect size
  power <- pwr.t.test(d = effect_size, n = DZsample_size, sig.level = alpha, power = NULL, type = c("one.sample"), alternative = c("two.sided"))
  
  # Append the result to the data frame
  DZresults <- rbind(DZresults, data.frame(EffectSize = effect_size, Power = power$power))
}

# Create the plot
DZ_pw1 <- ggplot(DZresults, aes(x = EffectSize, y = Power)) +
  geom_line() +
  geom_point() +
  labs(x = "Effect Size", y = "Power") +
  scale_x_continuous(breaks = seq(0, 0.5, by = 0.05)) +
  scale_y_continuous(breaks = seq(0, 1, by = 0.05)) +
  ggtitle("Sensitivity power analysis (N=2030)") +
  theme_bw() +
  theme(plot.title = element_text(size = 10))



# DZ twin differences sample size
DZdiffsample_size <- 1584 # Sample size of DZ twin differences

# Initialize an empty data frame to store results
DZdiff_results <- data.frame(EffectSize = numeric(0), Power = numeric(0))

# Loop through different effect sizes
for (i in 1:length(effect_sizes)) {
  effect_size <- effect_sizes[i]
  
  # Calculate power for the current effect size
  power <- pwr.t.test(d = effect_size, n = DZdiffsample_size, sig.level = alpha, power = NULL, type = c("one.sample"), alternative = c("two.sided"))
  
  # Append the result to the data frame
  DZdiff_results <- rbind(DZdiff_results, data.frame(EffectSize = effect_size, Power = power$power))
}

# Create the plot using ggplot2
DZ_pw2 <- ggplot(DZdiff_results, aes(x = EffectSize, y = Power)) +
  geom_line() +
  geom_point() +
  labs(x = "Effect Size", y = "Power") +
  scale_x_continuous(breaks = seq(0, 0.5, by = 0.05)) +
  scale_y_continuous(breaks = seq(0, 1, by = 0.05)) +
  ggtitle("Sensitivity power analysis (N=1584)") +
  theme_bw() +
  theme(plot.title = element_text(size = 10))

grid.arrange(DZ_pw1, DZ_pw2, ncol = 2)

```

\clearpage

## Full regression results

\noindent \textbf{The effects of EA PGI and education on immigration attitudes (corresponds to Figure 4 in the main text)} \setstretch{1.2}

```{=tex}
\begin{table}[H]
\centering
\begin{threeparttable}
\caption{The effects of EA PGI and education on immigration attitudes (DZ twins)}
\label{tab:DZ-PGI}
\scriptsize
\setlength{\tabcolsep}{1pt}

\begin{tabular}{lccccccc} \hline
 & (1) & (2) & (3) & (4) & (5) & (6) & (7) \\
VARIABLES & Avg index & Avg index & Avg index & IRT index & IRT index & IRT index & Education \\ \hline
 &  &  &  &  &  &  &  \\
Years of Education & 0.087*** &  & 0.083*** & 0.084*** &  & 0.079*** &  \\
 & (0.011) &  & (0.011) & (0.011) &  & (0.011) &  \\
EA PGI &  & 0.133*** & 0.095** &  & 0.139*** & 0.103** & 0.455*** \\
 &  & (0.034) & (0.033) &  & (0.034) & (0.034) & (0.071) \\
Constant & -1.021*** & 0.009 & -0.964*** & -0.989*** & 0.005 & -0.926*** & 11.795*** \\
 & (0.133) & (0.024) & (0.135) & (0.136) & (0.024) & (0.138) & (0.056) \\
 &  &  &  &  &  &  &  \\
Observations & 3,168 & 3,168 & 3,168 & 3,168 & 3,168 & 3,168 & 3,168 \\
R-squared & 0.631 & 0.620 & 0.633 & 0.626 & 0.616 & 0.628 & 0.729 \\
Sex & YES & YES & YES & YES & YES & YES & YES \\
 Twin FE & YES & YES & YES & YES & YES & YES & YES \\ \hline
 
\end{tabular}
Notes: Standard errors are clustered at twin-pair level. *** p$<$0.001, ** p$<$0.01, * p$<$0.05.
\end{threeparttable}
\end{table}
```

\noindent \textbf{The effects of EA PGI and education by attitude items}

```{=tex}
\begin{table}[H]
\centering
\begin{threeparttable}
\caption{The effects of EA PGI and education by attitude items (DZ twins)}
\label{tab:DZ-PGI-item}
\scriptsize
\setlength{\tabcolsep}{1pt}

\begin{tabular}{lcccccccccccc} \hline
 & (1) & (2) & (3) & (4) & (5) & (6) & (7) & (8) & (9) & (10) & (11) & (12) \\
VARIABLES & refugee & refugee & refugee & labor & labor & labor & culture & culture & culture & language & language & language \\ \hline
 &  &  &  &  &  &  &  &  &  &  &  &  \\
Years of Education & 0.088*** &  & 0.083*** & 0.078*** &  & 0.073*** & 0.044*** &  & 0.045*** & 0.078*** &  & 0.070*** \\
 & (0.014) &  & (0.014) & (0.012) &  & (0.012) & (0.012) &  & (0.012) & (0.015) &  & (0.015) \\
EA PGI &  & 0.134** & 0.096* &  & 0.127*** & 0.094** &  & 0.001 & -0.019 &  & 0.174*** & 0.142*** \\
 &  & (0.042) & (0.041) &  & (0.036) & (0.036) &  & (0.035) & (0.035) &  & (0.041) & (0.041) \\
Constant & 1.748*** & 2.784*** & 1.806*** & 1.705*** & 2.625*** & 1.762*** & 1.654*** & 2.170*** & 1.642*** & 1.512*** & 2.428*** & 1.597*** \\
 & (0.171) & (0.030) & (0.173) & (0.144) & (0.026) & (0.146) & (0.137) & (0.025) & (0.138) & (0.174) & (0.031) & (0.176) \\
 &  &  &  &  &  &  &  &  &  &  &  &  \\
Observations & 3,168 & 3,168 & 3,168 & 3,168 & 3,168 & 3,168 & 3,168 & 3,168 & 3,168 & 3,168 & 3,168 & 3,168 \\
R-squared & 0.608 & 0.601 & 0.610 & 0.585 & 0.577 & 0.587 & 0.577 & 0.573 & 0.577 & 0.575 & 0.572 & 0.578 \\
Sex & YES & YES & YES & YES & YES & YES & YES & YES & YES & YES & YES & YES \\
 Twin FE & YES & YES & YES & YES & YES & YES & YES & YES & YES & YES & YES & YES \\ \hline
\end{tabular}

Notes: Standard errors are clustered at twin-pair level. Dependent variables are on the original 1-5 scale. *** p$<$0.001, ** p$<$0.01, * p$<$0.05.
\end{threeparttable}
\end{table}
```

\clearpage

\noindent \textbf{Between family analysis of EA PGI, education, and immigration attitudes}

```{=tex}
\begin{table}[H]
\centering
\begin{threeparttable}
\caption{The effects of EA PGI and education on immigration attitudes (Between family)}
\label{tab:DZ-PGI}
\scriptsize
\setlength{\tabcolsep}{1pt}

\begin{tabular}{lccccccc} \hline
 & (1) & (2) & (3) & (4) & (5) & (6) & (7) \\
VARIABLES & Avg index & Avg index & Avg index & IRT index & IRT index & IRT index & Education \\ \hline
 &  &  &  &  &  &  &  \\
Years of Education & 0.109*** &  & 0.098*** & 0.115*** &  & 0.103*** &  \\
 & (0.007) &  & (0.007) & (0.007) &  & (0.007) &  \\
EA PGI &  & 0.170*** & 0.095*** &  & 0.188*** & 0.109*** & 0.766*** \\
 &  & (0.018) & (0.019) &  & (0.018) & (0.018) & (0.048) \\
Constant & -1.321*** & -0.064 & -1.193*** & -1.371*** & -0.038 & -1.224*** & 11.488*** \\
 & (0.128) & (0.099) & (0.130) & (0.130) & (0.101) & (0.132) & (0.275) \\
 &  &  &  &  &  &  &  \\
Observations & 3,168 & 3,168 & 3,168 & 3,168 & 3,168 & 3,168 & 3,168 \\
R-squared & 0.109 & 0.057 & 0.117 & 0.118 & 0.063 & 0.129 & 0.119 \\
Controls & YES & YES & YES & YES & YES & YES & YES \\
 Twin FE & NO & NO & NO & NO & NO & NO & NO \\ \hline
\end{tabular}
Notes: All models include controls for sex, birth-year dummies, interaction between sex and the birth-year dummies, the first twenty principal components of the genetic relationship matrix, and batch fixed effects. Standard errors are clustered at twin-pair level. *** p$<$0.001, ** p$<$0.01, * p$<$0.05.
\end{threeparttable}
\end{table}
```
\clearpage



## A simple test on mediators

To explore other potential mediators between EA PGI and immigration attitudes, I conducted some simple test with the following criteria: EA PGI should have significant effect on the mediator, and the mediator should plausibly affect immigration attitudes beyond education.The results reported here should be interpreted as only suggestive.

I focus on a psychological trait - locus of control (LOC) from available data. The LOC classifies individuals along a single dimension capturing the degree to which they feel like they control the outcome of events. Those with an internal locus of control believe they have command over their own fate, attributing achieved outcomes to their own hard work and abilities. Conversely, individuals with an external locus of control perceive outcomes as beyond their control. Individual LOC has been shown to predict immigration sentiments - those who feel in control are less hostile to immigrants [@harell_locus_2017]. There has also been extensive literature showing the connection between LOC and academic achievements (e.g., @findley_locus_1983). The LOC used is based on a 12-item forced-choice questionnaire [@rotter_psychological_1966], and the responses to each choice set are summed together. LOC is reversely coded so that higher scores are associated with higher internal locus of control. The variable is standardised.

The first column in the table below shows that EA PGI has significant effects on LOC. Column 2 confirms that higher personal control has positive effects on pro-immigration attitudes, while column 5 shows that LOC slightly attenuates the effects of EA PGI by 6%. Including education decrease the estimates for both EA PGI and LOC, but LOC remains a significant predictor for immigration attitudes. Comparing column 4 and 6 where education has been controlled for, the inclusion of LOC attenuates the effects of EA PGI by about 3%. The patterns for the IRT index are largely identical to those for the average index.

```{=tex}
\begin{table}[H]
\centering
\begin{threeparttable}
\caption{The effects of EA PGI and Locus of control on immigration attitudes (DZ twins)}
\label{tab:DZ-PGI-LOC}
\scriptsize
\setlength{\tabcolsep}{1pt}

\begin{tabular}{lccccccccccc} \hline
 & (1) & (2) & (3) & (4) & (5) & (6) & (7) & (8) & (9) & (10) & (11) \\
VARIABLES & LOC & Avg index & Avg index & Avg index & Avg index & Avg index & IRT index & IRT index & IRT index & IRT index & IRT index \\ \hline
 &  &  &  &  &  &  &  &  &  &  &  \\
EA PGI & 0.076* &  & 0.136** & 0.095* & 0.128** & 0.092* &  & 0.139*** & 0.100* & 0.132** & 0.097* \\
 & (0.038) &  & (0.041) & (0.041) & (0.041) & (0.040) &  & (0.042) & (0.041) & (0.042) & (0.041) \\
LOC &  & 0.100** &  &  & 0.094** & 0.075* & 0.096** &  &  & 0.090** & 0.072* \\
 &  & (0.031) &  &  & (0.030) & (0.030) & (0.030) &  &  & (0.030) & (0.030) \\
Years of Education &  &  &  & 0.084*** &  & 0.080*** &  &  & 0.082*** &  & 0.078*** \\
 &  &  &  & (0.014) &  & (0.014) &  &  & (0.014) &  & (0.014) \\
Constant & 0.090** & -0.011 & 0.001 & -1.002*** & -0.008 & -0.963*** & -0.019 & -0.008 & -0.979*** & -0.016 & -0.941*** \\
 & (0.028) & (0.027) & (0.027) & (0.166) & (0.027) & (0.164) & (0.028) & (0.028) & (0.168) & (0.027) & (0.167) \\
 &  &  &  &  &  &  &  &  &  &  &  \\
Observations & 2,296 & 2,296 & 2,296 & 2,296 & 2,296 & 2,296 & 2,296 & 2,296 & 2,296 & 2,296 & 2,296 \\
R-squared & 0.575 & 0.611 & 0.611 & 0.624 & 0.615 & 0.627 & 0.609 & 0.610 & 0.622 & 0.613 & 0.624 \\
Sex & YES & YES & YES & YES & YES & YES & YES & YES & YES & YES & YES \\
 Twin FE & YES & YES & YES & YES & YES & YES & YES & YES & YES & YES & YES \\ \hline
 
\end{tabular}
Notes: Standard errors are clustered at twin-pair level. LOC is a standardised variable based on a 12-item battery.  *** p$<$0.001, ** p$<$0.01, * p$<$0.05.
\end{threeparttable}
\end{table}
```


\noindent \textbf{Locus of Control items} \vspace{-2pt} \singlespacing

```{=tex}
\noindent \textit{For each pair of statement choose the one that best describes how you feel.}

\footnotesize
1. Many of the unhappy things in people's lives are partly due to bad luck.

2. People's misfortunes result from the mistakes they make.
\vspace{6pt}

1. An important reason why there are wars is that people are not interested enough in politics.

2. There will always be war no matter what many efforts people make to prevent it.
\vspace{6pt}

1. In the long run, people get the respect they deserve in this world.

2. Unfortunately, an individual's worth often passes unrecognized no matter how hard he tries.
\vspace{6pt}

1. The idea that teachers are unfair to students is nonsense.

2. Most students don't realize the extent to which their grades are affected by accidental happenings.
\vspace{6pt}

1. Without the right luck, one cannot be an effective leader.

2. Capable people who fail to became leaders have not taken advantage of their opportunities.
\vspace{6pt}

1. No matter how hard you try, some people just don't like you.

2. People who can't get others to like them don't understand how to get along with others.
\vspace{6pt}

1. What happens, happens.

2. Trusting fate has never worked as well for me as deciding on a particular plan of action.
\vspace{6pt}

1. Becoming a success is a matter of hard work; luck has little or nothing to do with.

2. Getting a good job depends mainly on being in the right place at the right time.
\vspace{6pt}

1. An average citizen can influence government decisions.

2. The world is ruled by a few people in power and there is not much that the little man can do about.
\vspace{6pt}

1. When I make plans, I am almost certain that I can make them work.

2. It is not always wise to plan too far ahead because many things turn out to be a matter of luck anyway.
\vspace{6pt}

1. In my case, getting what I want has little or nothing to do with luck.

2. Many times we might just as well decide what to do by flipping a coin.
\vspace{6pt}

1. What happens to me is my own doing. 

2. Sometimes I feel that I don't have enough control over the direction my life is taking.
```



# References

\small
